File indexing completed on 2023-01-12 08:55:39
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051 #include "SimG4Core/CustomPhysics/interface/FullModelReactionDynamics.h"
0052 #include "G4AntiProton.hh"
0053 #include "G4AntiNeutron.hh"
0054 #include "Randomize.hh"
0055 #include <iostream>
0056 #include "G4ParticleTable.hh"
0057
0058 using namespace CLHEP;
0059
0060 G4bool FullModelReactionDynamics::GenerateXandPt(
0061 G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
0062 G4int &vecLen,
0063 G4ReactionProduct &modifiedOriginal,
0064 const G4HadProjectile *originalIncident,
0065 G4ReactionProduct ¤tParticle,
0066 G4ReactionProduct &targetParticle,
0067 const G4Nucleus &targetNucleus,
0068 G4bool &incidentHasChanged,
0069 G4bool &targetHasChanged,
0070 G4bool leadFlag,
0071 G4ReactionProduct &leadingStrangeParticle) {
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085 if (vecLen == 0)
0086 return false;
0087
0088 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
0089
0090
0091 G4ParticleDefinition *aProton = G4Proton::Proton();
0092 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
0093 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
0094 G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
0095 G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
0096 G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
0097 G4ParticleDefinition *aKaonZeroS = G4KaonZeroShort::KaonZeroShort();
0098 G4ParticleDefinition *aKaonZeroL = G4KaonZeroLong::KaonZeroLong();
0099
0100 G4int i, l;
0101
0102 G4bool veryForward = false;
0103
0104 const G4double ekOriginal = modifiedOriginal.GetKineticEnergy() / GeV;
0105 const G4double etOriginal = modifiedOriginal.GetTotalEnergy() / GeV;
0106 const G4double mOriginal = modifiedOriginal.GetMass() / GeV;
0107 const G4double pOriginal = modifiedOriginal.GetMomentum().mag() / GeV;
0108 G4double targetMass = targetParticle.GetDefinition()->GetPDGMass() / GeV;
0109 G4double centerofmassEnergy =
0110 std::sqrt(mOriginal * mOriginal + targetMass * targetMass + 2.0 * targetMass * etOriginal);
0111 G4double currentMass = currentParticle.GetMass() / GeV;
0112 targetMass = targetParticle.GetMass() / GeV;
0113
0114
0115
0116
0117 for (i = 0; i < vecLen; ++i) {
0118 G4int itemp = G4int(G4UniformRand() * vecLen);
0119 G4ReactionProduct pTemp = *vec[itemp];
0120 *vec[itemp] = *vec[i];
0121 *vec[i] = pTemp;
0122
0123 }
0124
0125 if (currentMass == 0.0 && targetMass == 0.0)
0126 {
0127
0128 G4double ek = currentParticle.GetKineticEnergy();
0129 G4ThreeVector m = currentParticle.GetMomentum();
0130 currentParticle = *vec[0];
0131 targetParticle = *vec[1];
0132 for (i = 0; i < (vecLen - 2); ++i)
0133 *vec[i] = *vec[i + 2];
0134 G4ReactionProduct *temp = vec[vecLen - 1];
0135 delete temp;
0136 temp = vec[vecLen - 2];
0137 delete temp;
0138 vecLen -= 2;
0139 currentMass = currentParticle.GetMass() / GeV;
0140 targetMass = targetParticle.GetMass() / GeV;
0141 incidentHasChanged = true;
0142 targetHasChanged = true;
0143 currentParticle.SetKineticEnergy(ek);
0144 currentParticle.SetMomentum(m);
0145 veryForward = true;
0146 }
0147 const G4double atomicWeight = targetNucleus.GetN_asInt();
0148 const G4double atomicNumber = targetNucleus.GetZ_asInt();
0149 const G4double protonMass = aProton->GetPDGMass() / MeV;
0150 if ((originalIncident->GetDefinition() == aKaonMinus || originalIncident->GetDefinition() == aKaonZeroL ||
0151 originalIncident->GetDefinition() == aKaonZeroS || originalIncident->GetDefinition() == aKaonPlus) &&
0152 G4UniformRand() >= 0.7) {
0153 G4ReactionProduct temp = currentParticle;
0154 currentParticle = targetParticle;
0155 targetParticle = temp;
0156 incidentHasChanged = true;
0157 targetHasChanged = true;
0158 currentMass = currentParticle.GetMass() / GeV;
0159 targetMass = targetParticle.GetMass() / GeV;
0160 }
0161 const G4double afc = std::min(0.75,
0162 0.312 + 0.200 * std::log(std::log(centerofmassEnergy * centerofmassEnergy)) +
0163 std::pow(centerofmassEnergy * centerofmassEnergy, 1.5) / 6000.0);
0164
0165
0166
0167 G4double freeEnergy = centerofmassEnergy - currentMass - targetMass;
0168
0169 if (freeEnergy < 0) {
0170 G4cout << "Free energy < 0!" << G4endl;
0171 G4cout << "E_CMS = " << centerofmassEnergy << " GeV" << G4endl;
0172 G4cout << "m_curr = " << currentMass << " GeV" << G4endl;
0173 G4cout << "m_orig = " << mOriginal << " GeV" << G4endl;
0174 G4cout << "m_targ = " << targetMass << " GeV" << G4endl;
0175 G4cout << "E_free = " << freeEnergy << " GeV" << G4endl;
0176 }
0177
0178 G4double forwardEnergy = freeEnergy / 2.;
0179 G4int forwardCount = 1;
0180
0181 G4double backwardEnergy = freeEnergy / 2.;
0182 G4int backwardCount = 1;
0183 if (veryForward) {
0184 if (currentParticle.GetSide() == -1) {
0185 forwardEnergy += currentMass;
0186 forwardCount--;
0187 backwardEnergy -= currentMass;
0188 backwardCount++;
0189 }
0190 if (targetParticle.GetSide() != -1) {
0191 backwardEnergy += targetMass;
0192 backwardCount--;
0193 forwardEnergy -= targetMass;
0194 forwardCount++;
0195 }
0196 }
0197 for (i = 0; i < vecLen; ++i) {
0198 if (vec[i]->GetSide() == -1) {
0199 ++backwardCount;
0200 backwardEnergy -= vec[i]->GetMass() / GeV;
0201 } else {
0202 ++forwardCount;
0203 forwardEnergy -= vec[i]->GetMass() / GeV;
0204 }
0205 }
0206
0207
0208
0209
0210
0211 G4double xtarg;
0212 if (centerofmassEnergy < (2.0 + G4UniformRand()))
0213 xtarg = afc * (std::pow(atomicWeight, 0.33) - 1.0) * (2.0 * backwardCount + vecLen + 2) / 2.0;
0214 else
0215 xtarg = afc * (std::pow(atomicWeight, 0.33) - 1.0) * (2.0 * backwardCount);
0216 if (xtarg <= 0.0)
0217 xtarg = 0.01;
0218 G4int nuclearExcitationCount = Poisson(xtarg);
0219 if (atomicWeight < 1.0001)
0220 nuclearExcitationCount = 0;
0221 G4int extraNucleonCount = 0;
0222
0223 if (nuclearExcitationCount > 0) {
0224 const G4double nucsup[] = {1.00, 0.7, 0.5, 0.4, 0.35, 0.3};
0225 const G4double psup[] = {3., 6., 20., 50., 100., 1000.};
0226 G4int momentumBin = 0;
0227 while ((momentumBin < 6) && (modifiedOriginal.GetTotalMomentum() / GeV > psup[momentumBin]))
0228 ++momentumBin;
0229 momentumBin = std::min(5, momentumBin);
0230
0231
0232
0233
0234 for (i = 0; i < nuclearExcitationCount; ++i) {
0235 G4ReactionProduct *pVec = new G4ReactionProduct();
0236 if (G4UniformRand() < nucsup[momentumBin]) {
0237 if (G4UniformRand() > 1.0 - atomicNumber / atomicWeight)
0238 pVec->SetDefinition(aProton);
0239 else
0240 pVec->SetDefinition(aNeutron);
0241 pVec->SetSide(-2);
0242 ++extraNucleonCount;
0243 backwardEnergy += pVec->GetMass() / GeV;
0244
0245 } else {
0246 G4double ran = G4UniformRand();
0247 if (ran < 0.3181)
0248 pVec->SetDefinition(aPiPlus);
0249 else if (ran < 0.6819)
0250 pVec->SetDefinition(aPiZero);
0251 else
0252 pVec->SetDefinition(aPiMinus);
0253 pVec->SetSide(-1);
0254 }
0255 pVec->SetNewlyAdded(true);
0256 vec.SetElement(vecLen++, pVec);
0257
0258 backwardEnergy -= pVec->GetMass() / GeV;
0259 ++backwardCount;
0260 }
0261 }
0262
0263
0264
0265 G4int is, iskip;
0266 while (forwardEnergy <= 0.0)
0267 {
0268
0269 iskip = G4int(G4UniformRand() * forwardCount) + 1;
0270 is = 0;
0271 G4int forwardParticlesLeft = 0;
0272 for (i = (vecLen - 1); i >= 0; --i) {
0273 if (vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled()) {
0274 forwardParticlesLeft = 1;
0275 if (++is == iskip) {
0276 forwardEnergy += vec[i]->GetMass() / GeV;
0277 for (G4int j = i; j < (vecLen - 1); j++)
0278 *vec[j] = *vec[j + 1];
0279 --forwardCount;
0280 G4ReactionProduct *temp = vec[vecLen - 1];
0281 delete temp;
0282 if (--vecLen == 0)
0283 return false;
0284 break;
0285 }
0286 }
0287 }
0288
0289 if (forwardParticlesLeft == 0) {
0290 forwardEnergy += currentParticle.GetMass() / GeV;
0291 currentParticle.SetDefinitionAndUpdateE(targetParticle.GetDefinition());
0292 targetParticle.SetDefinitionAndUpdateE(vec[0]->GetDefinition());
0293
0294 --forwardCount;
0295 for (G4int j = 0; j < (vecLen - 1); ++j)
0296 *vec[j] = *vec[j + 1];
0297 G4ReactionProduct *temp = vec[vecLen - 1];
0298 delete temp;
0299 if (--vecLen == 0)
0300 return false;
0301 break;
0302 }
0303 }
0304
0305 while (backwardEnergy <= 0.0)
0306 {
0307
0308 iskip = G4int(G4UniformRand() * backwardCount) + 1;
0309 is = 0;
0310 G4int backwardParticlesLeft = 0;
0311 for (i = (vecLen - 1); i >= 0; --i) {
0312 if (vec[i]->GetSide() < 0 && vec[i]->GetMayBeKilled()) {
0313 backwardParticlesLeft = 1;
0314 if (++is == iskip)
0315 {
0316 if (vec[i]->GetSide() == -2) {
0317 --extraNucleonCount;
0318
0319 backwardEnergy -= vec[i]->GetTotalEnergy() / GeV;
0320 }
0321 backwardEnergy += vec[i]->GetTotalEnergy() / GeV;
0322 for (G4int j = i; j < (vecLen - 1); ++j)
0323 *vec[j] = *vec[j + 1];
0324 --backwardCount;
0325 G4ReactionProduct *temp = vec[vecLen - 1];
0326 delete temp;
0327 if (--vecLen == 0)
0328 return false;
0329 break;
0330 }
0331 }
0332 }
0333
0334 if (backwardParticlesLeft == 0) {
0335 backwardEnergy += targetParticle.GetMass() / GeV;
0336 targetParticle = *vec[0];
0337 --backwardCount;
0338 for (G4int j = 0; j < (vecLen - 1); ++j)
0339 *vec[j] = *vec[j + 1];
0340 G4ReactionProduct *temp = vec[vecLen - 1];
0341 delete temp;
0342 if (--vecLen == 0)
0343 return false;
0344 break;
0345 }
0346 }
0347
0348
0349
0350
0351
0352 G4ReactionProduct pseudoParticle[10];
0353 for (i = 0; i < 10; ++i)
0354 pseudoParticle[i].SetZero();
0355
0356 pseudoParticle[0].SetMass(mOriginal * GeV);
0357 pseudoParticle[0].SetMomentum(0.0, 0.0, pOriginal * GeV);
0358 pseudoParticle[0].SetTotalEnergy(std::sqrt(pOriginal * pOriginal + mOriginal * mOriginal) * GeV);
0359
0360 pseudoParticle[1].SetMass(protonMass * MeV);
0361 pseudoParticle[1].SetTotalEnergy(protonMass * MeV);
0362
0363 pseudoParticle[3].SetMass(protonMass * (1 + extraNucleonCount) * MeV);
0364 pseudoParticle[3].SetTotalEnergy(protonMass * (1 + extraNucleonCount) * MeV);
0365
0366 pseudoParticle[8].SetMomentum(1.0 * GeV, 0.0, 0.0);
0367
0368 pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
0369 pseudoParticle[3] = pseudoParticle[3] + pseudoParticle[0];
0370
0371 pseudoParticle[0].Lorentz(pseudoParticle[0], pseudoParticle[2]);
0372 pseudoParticle[1].Lorentz(pseudoParticle[1], pseudoParticle[2]);
0373
0374 G4double dndl[20];
0375
0376
0377
0378
0379 G4double aspar, pt, et, x, pp, pp1, rthnve, phinve, rmb, wgt;
0380 G4int innerCounter, outerCounter;
0381 G4bool eliminateThisParticle, resetEnergies, constantCrossSection;
0382
0383 G4double forwardKinetic = 0.0, backwardKinetic = 0.0;
0384
0385
0386
0387
0388
0389 G4double binl[20] = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
0390 1.0, 1.11, 1.25, 1.43, 1.67, 2.0, 2.5, 3.33, 5.00, 10.00};
0391 G4int backwardNucleonCount = 0;
0392 G4double totalEnergy, kineticEnergy, vecMass;
0393
0394 for (i = (vecLen - 1); i >= 0; --i) {
0395 G4double phi = G4UniformRand() * twopi;
0396 if (vec[i]->GetNewlyAdded())
0397 {
0398 if (vec[i]->GetSide() == -2)
0399 {
0400 if (backwardNucleonCount < 18) {
0401 if (vec[i]->GetDefinition() == G4PionMinus::PionMinus() ||
0402 vec[i]->GetDefinition() == G4PionPlus::PionPlus() || vec[i]->GetDefinition() == G4PionZero::PionZero()) {
0403 for (G4int i = 0; i < vecLen; i++)
0404 delete vec[i];
0405 vecLen = 0;
0406 G4ExceptionDescription ed;
0407 ed << "FullModelReactionDynamics::GenerateXandPt : a pion has been counted as a backward nucleon";
0408 G4Exception("FullModelReactionDynamics::GenerateXandPt", "had064", FatalException, ed);
0409 }
0410 vec[i]->SetSide(-3);
0411 ++backwardNucleonCount;
0412 continue;
0413 }
0414 }
0415 }
0416
0417
0418
0419
0420 vecMass = vec[i]->GetMass() / GeV;
0421 G4double ran = -std::log(1.0 - G4UniformRand()) / 3.5;
0422 if (vec[i]->GetSide() == -2)
0423 {
0424 if (vec[i]->GetDefinition() == aKaonMinus || vec[i]->GetDefinition() == aKaonZeroL ||
0425 vec[i]->GetDefinition() == aKaonZeroS || vec[i]->GetDefinition() == aKaonPlus ||
0426 vec[i]->GetDefinition() == aPiMinus || vec[i]->GetDefinition() == aPiZero ||
0427 vec[i]->GetDefinition() == aPiPlus) {
0428 aspar = 0.75;
0429 pt = std::sqrt(std::pow(ran, 1.7));
0430 } else {
0431 aspar = 0.20;
0432 pt = std::sqrt(std::pow(ran, 1.2));
0433 }
0434 } else {
0435 if (vec[i]->GetDefinition() == aPiMinus || vec[i]->GetDefinition() == aPiZero ||
0436 vec[i]->GetDefinition() == aPiPlus) {
0437 aspar = 0.75;
0438 pt = std::sqrt(std::pow(ran, 1.7));
0439 } else if (vec[i]->GetDefinition() == aKaonMinus || vec[i]->GetDefinition() == aKaonZeroL ||
0440 vec[i]->GetDefinition() == aKaonZeroS || vec[i]->GetDefinition() == aKaonPlus) {
0441 aspar = 0.70;
0442 pt = std::sqrt(std::pow(ran, 1.7));
0443 } else {
0444 aspar = 0.65;
0445 pt = std::sqrt(std::pow(ran, 1.5));
0446 }
0447 }
0448 pt = std::max(0.001, pt);
0449 vec[i]->SetMomentum(pt * std::cos(phi) * GeV, pt * std::sin(phi) * GeV);
0450 for (G4int j = 0; j < 20; ++j)
0451 binl[j] = j / (19. * pt);
0452 if (vec[i]->GetSide() > 0)
0453 et = pseudoParticle[0].GetTotalEnergy() / GeV;
0454 else
0455 et = pseudoParticle[1].GetTotalEnergy() / GeV;
0456 dndl[0] = 0.0;
0457
0458
0459
0460 outerCounter = 0;
0461 eliminateThisParticle = true;
0462 resetEnergies = true;
0463 while (++outerCounter < 3) {
0464 for (l = 1; l < 20; ++l) {
0465 x = (binl[l] + binl[l - 1]) / 2.;
0466 pt = std::max(0.001, pt);
0467 if (x > 1.0 / pt)
0468 dndl[l] += dndl[l - 1];
0469 else
0470 dndl[l] = et * aspar / std::sqrt(std::pow((1. + aspar * x * aspar * x), 3)) * (binl[l] - binl[l - 1]) /
0471 std::sqrt(pt * x * et * pt * x * et + pt * pt + vecMass * vecMass) +
0472 dndl[l - 1];
0473 }
0474 innerCounter = 0;
0475 vec[i]->SetMomentum(pt * std::cos(phi) * GeV, pt * std::sin(phi) * GeV);
0476
0477
0478
0479 while (++innerCounter < 7) {
0480 ran = G4UniformRand() * dndl[19];
0481 l = 1;
0482 while ((ran >= dndl[l]) && (l < 20))
0483 l++;
0484 l = std::min(19, l);
0485 x = std::min(1.0, pt * (binl[l - 1] + G4UniformRand() * (binl[l] - binl[l - 1]) / 2.));
0486 if (vec[i]->GetSide() < 0)
0487 x *= -1.;
0488 vec[i]->SetMomentum(x * et * GeV);
0489 totalEnergy = std::sqrt(x * et * x * et + pt * pt + vecMass * vecMass);
0490 vec[i]->SetTotalEnergy(totalEnergy * GeV);
0491 kineticEnergy = vec[i]->GetKineticEnergy() / GeV;
0492 if (vec[i]->GetSide() > 0)
0493 {
0494 if ((forwardKinetic + kineticEnergy) < 0.95 * forwardEnergy) {
0495 pseudoParticle[4] = pseudoParticle[4] + (*vec[i]);
0496 forwardKinetic += kineticEnergy;
0497 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
0498 pseudoParticle[6].SetMomentum(0.0);
0499 phi = pseudoParticle[6].Angle(pseudoParticle[8]);
0500 if (pseudoParticle[6].GetMomentum().y() / MeV < 0.0)
0501 phi = twopi - phi;
0502 phi += pi + normal() * pi / 12.0;
0503 if (phi > twopi)
0504 phi -= twopi;
0505 if (phi < 0.0)
0506 phi = twopi - phi;
0507 outerCounter = 2;
0508 eliminateThisParticle = false;
0509 resetEnergies = false;
0510 break;
0511 }
0512 if (innerCounter > 5)
0513 break;
0514 if (backwardEnergy >= vecMass)
0515 {
0516 vec[i]->SetSide(-1);
0517 forwardEnergy += vecMass;
0518 backwardEnergy -= vecMass;
0519 ++backwardCount;
0520 }
0521 } else {
0522 G4double xxx = 0.95 + 0.05 * extraNucleonCount / 20.0;
0523 if ((backwardKinetic + kineticEnergy) < xxx * backwardEnergy) {
0524 pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
0525 backwardKinetic += kineticEnergy;
0526 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
0527 pseudoParticle[6].SetMomentum(0.0);
0528 phi = pseudoParticle[6].Angle(pseudoParticle[8]);
0529 if (pseudoParticle[6].GetMomentum().y() / MeV < 0.0)
0530 phi = twopi - phi;
0531 phi += pi + normal() * pi / 12.0;
0532 if (phi > twopi)
0533 phi -= twopi;
0534 if (phi < 0.0)
0535 phi = twopi - phi;
0536 outerCounter = 2;
0537 eliminateThisParticle = false;
0538 resetEnergies = false;
0539 break;
0540 }
0541 if (innerCounter > 5)
0542 break;
0543 if (forwardEnergy >= vecMass)
0544 {
0545 vec[i]->SetSide(1);
0546 forwardEnergy -= vecMass;
0547 backwardEnergy += vecMass;
0548 backwardCount--;
0549 }
0550 }
0551 G4ThreeVector momentum = vec[i]->GetMomentum();
0552 vec[i]->SetMomentum(momentum.x() * 0.9, momentum.y() * 0.9);
0553 pt *= 0.9;
0554 dndl[19] *= 0.9;
0555 }
0556 if (resetEnergies) {
0557
0558
0559
0560
0561
0562 forwardKinetic = 0.0;
0563 backwardKinetic = 0.0;
0564 pseudoParticle[4].SetZero();
0565 pseudoParticle[5].SetZero();
0566 for (l = i + 1; l < vecLen; ++l) {
0567 if (vec[l]->GetSide() > 0 || vec[l]->GetDefinition() == aKaonMinus || vec[l]->GetDefinition() == aKaonZeroL ||
0568 vec[l]->GetDefinition() == aKaonZeroS || vec[l]->GetDefinition() == aKaonPlus ||
0569 vec[l]->GetDefinition() == aPiMinus || vec[l]->GetDefinition() == aPiZero ||
0570 vec[l]->GetDefinition() == aPiPlus) {
0571 G4double tempMass = vec[l]->GetMass() / MeV;
0572 totalEnergy = 0.95 * vec[l]->GetTotalEnergy() / MeV + 0.05 * tempMass;
0573 totalEnergy = std::max(tempMass, totalEnergy);
0574 vec[l]->SetTotalEnergy(totalEnergy * MeV);
0575 pp = std::sqrt(std::abs(totalEnergy * totalEnergy - tempMass * tempMass));
0576 pp1 = vec[l]->GetMomentum().mag() / MeV;
0577 if (pp1 < 1.0e-6 * GeV) {
0578 G4double rthnve = pi * G4UniformRand();
0579 G4double phinve = twopi * G4UniformRand();
0580 G4double srth = std::sin(rthnve);
0581 vec[l]->SetMomentum(
0582 pp * srth * std::cos(phinve) * MeV, pp * srth * std::sin(phinve) * MeV, pp * std::cos(rthnve) * MeV);
0583 } else {
0584 vec[l]->SetMomentum(vec[l]->GetMomentum() * (pp / pp1));
0585 }
0586 G4double px = vec[l]->GetMomentum().x() / MeV;
0587 G4double py = vec[l]->GetMomentum().y() / MeV;
0588 pt = std::max(1.0, std::sqrt(px * px + py * py)) / GeV;
0589 if (vec[l]->GetSide() > 0) {
0590 forwardKinetic += vec[l]->GetKineticEnergy() / GeV;
0591 pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
0592 } else {
0593 backwardKinetic += vec[l]->GetKineticEnergy() / GeV;
0594 pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
0595 }
0596 }
0597 }
0598 }
0599 }
0600
0601 if (eliminateThisParticle && vec[i]->GetMayBeKilled())
0602 {
0603 if (vec[i]->GetSide() > 0) {
0604 --forwardCount;
0605 forwardEnergy += vecMass;
0606 } else {
0607 if (vec[i]->GetSide() == -2) {
0608 --extraNucleonCount;
0609
0610 backwardEnergy -= vecMass;
0611 }
0612 --backwardCount;
0613 backwardEnergy += vecMass;
0614 }
0615 for (G4int j = i; j < (vecLen - 1); ++j)
0616 *vec[j] = *vec[j + 1];
0617 G4ReactionProduct *temp = vec[vecLen - 1];
0618 delete temp;
0619
0620 if (--vecLen == 0)
0621 return false;
0622 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
0623 pseudoParticle[6].SetMomentum(0.0);
0624 }
0625 }
0626
0627
0628
0629
0630
0631
0632 G4double phi = G4UniformRand() * twopi;
0633 G4double ran = -std::log(1.0 - G4UniformRand());
0634 if (currentParticle.GetDefinition() == aPiMinus || currentParticle.GetDefinition() == aPiZero ||
0635 currentParticle.GetDefinition() == aPiPlus) {
0636 aspar = 0.60;
0637 pt = std::sqrt(std::pow(ran / 6.0, 1.7));
0638 } else if (currentParticle.GetDefinition() == aKaonMinus || currentParticle.GetDefinition() == aKaonZeroL ||
0639 currentParticle.GetDefinition() == aKaonZeroS || currentParticle.GetDefinition() == aKaonPlus) {
0640 aspar = 0.50;
0641 pt = std::sqrt(std::pow(ran / 5.0, 1.4));
0642 } else {
0643 aspar = 0.40;
0644 pt = std::sqrt(std::pow(ran / 4.0, 1.2));
0645 }
0646 for (G4int j = 0; j < 20; ++j)
0647 binl[j] = j / (19. * pt);
0648 currentParticle.SetMomentum(pt * std::cos(phi) * GeV, pt * std::sin(phi) * GeV);
0649 et = pseudoParticle[0].GetTotalEnergy() / GeV;
0650 dndl[0] = 0.0;
0651 vecMass = currentParticle.GetMass() / GeV;
0652 for (l = 1; l < 20; ++l) {
0653 x = (binl[l] + binl[l - 1]) / 2.;
0654 if (x > 1.0 / pt)
0655 dndl[l] += dndl[l - 1];
0656 else
0657 dndl[l] = aspar / std::sqrt(std::pow((1. + sqr(aspar * x)), 3)) * (binl[l] - binl[l - 1]) * et /
0658 std::sqrt(pt * x * et * pt * x * et + pt * pt + vecMass * vecMass) +
0659 dndl[l - 1];
0660 }
0661 ran = G4UniformRand() * dndl[19];
0662 l = 1;
0663 while ((ran > dndl[l]) && (l < 20))
0664 l++;
0665 l = std::min(19, l);
0666 x = std::min(1.0, pt * (binl[l - 1] + G4UniformRand() * (binl[l] - binl[l - 1]) / 2.));
0667 currentParticle.SetMomentum(x * et * GeV);
0668 if (forwardEnergy < forwardKinetic)
0669 totalEnergy = vecMass + 0.04 * std::fabs(normal());
0670 else
0671 totalEnergy = vecMass + forwardEnergy - forwardKinetic;
0672 currentParticle.SetTotalEnergy(totalEnergy * GeV);
0673 pp = std::sqrt(std::abs(totalEnergy * totalEnergy - vecMass * vecMass)) * GeV;
0674 pp1 = currentParticle.GetMomentum().mag() / MeV;
0675 if (pp1 < 1.0e-6 * GeV) {
0676 G4double rthnve = pi * G4UniformRand();
0677 G4double phinve = twopi * G4UniformRand();
0678 G4double srth = std::sin(rthnve);
0679 currentParticle.SetMomentum(
0680 pp * srth * std::cos(phinve) * MeV, pp * srth * std::sin(phinve) * MeV, pp * std::cos(rthnve) * MeV);
0681 } else {
0682 currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
0683 }
0684 pseudoParticle[4] = pseudoParticle[4] + currentParticle;
0685
0686
0687
0688
0689 if (backwardNucleonCount < 18) {
0690 targetParticle.SetSide(-3);
0691 ++backwardNucleonCount;
0692 } else {
0693
0694
0695
0696 vecMass = targetParticle.GetMass() / GeV;
0697 ran = -std::log(1.0 - G4UniformRand());
0698 aspar = 0.40;
0699 pt = std::max(0.001, std::sqrt(std::pow(ran / 4.0, 1.2)));
0700 targetParticle.SetMomentum(pt * std::cos(phi) * GeV, pt * std::sin(phi) * GeV);
0701 for (G4int j = 0; j < 20; ++j)
0702 binl[j] = (j - 1.) / (19. * pt);
0703 et = pseudoParticle[1].GetTotalEnergy() / GeV;
0704 dndl[0] = 0.0;
0705 outerCounter = 0;
0706 resetEnergies = true;
0707 while (++outerCounter < 3)
0708 {
0709 for (l = 1; l < 20; ++l) {
0710 x = (binl[l] + binl[l - 1]) / 2.;
0711 if (x > 1.0 / pt)
0712 dndl[l] += dndl[l - 1];
0713 else
0714 dndl[l] = aspar / std::sqrt(std::pow((1. + aspar * x * aspar * x), 3)) * (binl[l] - binl[l - 1]) * et /
0715 std::sqrt(pt * x * et * pt * x * et + pt * pt + vecMass * vecMass) +
0716 dndl[l - 1];
0717 }
0718 innerCounter = 0;
0719 while (++innerCounter < 7)
0720 {
0721 l = 1;
0722 ran = G4UniformRand() * dndl[19];
0723 while ((ran >= dndl[l]) && (l < 20))
0724 l++;
0725 l = std::min(19, l);
0726 x = std::min(1.0, pt * (binl[l - 1] + G4UniformRand() * (binl[l] - binl[l - 1]) / 2.));
0727 if (targetParticle.GetSide() < 0)
0728 x *= -1.;
0729 targetParticle.SetMomentum(x * et * GeV);
0730 totalEnergy = std::sqrt(x * et * x * et + pt * pt + vecMass * vecMass);
0731 targetParticle.SetTotalEnergy(totalEnergy * GeV);
0732 if (targetParticle.GetSide() < 0) {
0733 G4double xxx = 0.95 + 0.05 * extraNucleonCount / 20.0;
0734 if ((backwardKinetic + totalEnergy - vecMass) < xxx * backwardEnergy) {
0735 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
0736 backwardKinetic += totalEnergy - vecMass;
0737 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
0738 pseudoParticle[6].SetMomentum(0.0);
0739 outerCounter = 2;
0740 resetEnergies = false;
0741 break;
0742 }
0743 if (innerCounter > 5)
0744 break;
0745 if (forwardEnergy >= vecMass)
0746 {
0747 targetParticle.SetSide(1);
0748 forwardEnergy -= vecMass;
0749 backwardEnergy += vecMass;
0750 --backwardCount;
0751 }
0752 G4ThreeVector momentum = targetParticle.GetMomentum();
0753 targetParticle.SetMomentum(momentum.x() * 0.9, momentum.y() * 0.9);
0754 pt *= 0.9;
0755 dndl[19] *= 0.9;
0756 } else
0757 {
0758 if (forwardEnergy < forwardKinetic)
0759 totalEnergy = vecMass + 0.04 * std::fabs(normal());
0760 else
0761 totalEnergy = vecMass + forwardEnergy - forwardKinetic;
0762 targetParticle.SetTotalEnergy(totalEnergy * GeV);
0763 pp = std::sqrt(std::abs(totalEnergy * totalEnergy - vecMass * vecMass)) * GeV;
0764 pp1 = targetParticle.GetMomentum().mag() / MeV;
0765 if (pp1 < 1.0e-6 * GeV) {
0766 G4double rthnve = pi * G4UniformRand();
0767 G4double phinve = twopi * G4UniformRand();
0768 G4double srth = std::sin(rthnve);
0769 targetParticle.SetMomentum(
0770 pp * srth * std::cos(phinve) * MeV, pp * srth * std::sin(phinve) * MeV, pp * std::cos(rthnve) * MeV);
0771 } else
0772 targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
0773
0774 pseudoParticle[4] = pseudoParticle[4] + targetParticle;
0775 outerCounter = 2;
0776
0777 resetEnergies = false;
0778 break;
0779 }
0780 }
0781 if (resetEnergies) {
0782
0783
0784
0785
0786
0787 forwardKinetic = backwardKinetic = 0.0;
0788 pseudoParticle[4].SetZero();
0789 pseudoParticle[5].SetZero();
0790 for (l = 0; l < vecLen; ++l)
0791 {
0792 if (vec[l]->GetSide() > 0 || vec[l]->GetDefinition() == aKaonMinus || vec[l]->GetDefinition() == aKaonZeroL ||
0793 vec[l]->GetDefinition() == aKaonZeroS || vec[l]->GetDefinition() == aKaonPlus ||
0794 vec[l]->GetDefinition() == aPiMinus || vec[l]->GetDefinition() == aPiZero ||
0795 vec[l]->GetDefinition() == aPiPlus) {
0796 G4double tempMass = vec[l]->GetMass() / GeV;
0797 totalEnergy = std::max(tempMass, 0.95 * vec[l]->GetTotalEnergy() / GeV + 0.05 * tempMass);
0798 vec[l]->SetTotalEnergy(totalEnergy * GeV);
0799 pp = std::sqrt(std::abs(totalEnergy * totalEnergy - tempMass * tempMass)) * GeV;
0800 pp1 = vec[l]->GetMomentum().mag() / MeV;
0801 if (pp1 < 1.0e-6 * GeV) {
0802 G4double rthnve = pi * G4UniformRand();
0803 G4double phinve = twopi * G4UniformRand();
0804 G4double srth = std::sin(rthnve);
0805 vec[l]->SetMomentum(
0806 pp * srth * std::cos(phinve) * MeV, pp * srth * std::sin(phinve) * MeV, pp * std::cos(rthnve) * MeV);
0807 } else
0808 vec[l]->SetMomentum(vec[l]->GetMomentum() * (pp / pp1));
0809
0810 pt = std::max(0.001 * GeV,
0811 std::sqrt(sqr(vec[l]->GetMomentum().x() / MeV) + sqr(vec[l]->GetMomentum().y() / MeV))) /
0812 GeV;
0813 if (vec[l]->GetSide() > 0) {
0814 forwardKinetic += vec[l]->GetKineticEnergy() / GeV;
0815 pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
0816 } else {
0817 backwardKinetic += vec[l]->GetKineticEnergy() / GeV;
0818 pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
0819 }
0820 }
0821 }
0822 }
0823 }
0824 }
0825
0826
0827
0828
0829 pseudoParticle[6].Lorentz(pseudoParticle[3], pseudoParticle[2]);
0830 pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[4];
0831 pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[5];
0832 if (backwardNucleonCount == 1)
0833 {
0834 G4double ekin = std::min(backwardEnergy - backwardKinetic, centerofmassEnergy / 2.0 - protonMass / GeV);
0835 if (ekin < 0.04)
0836 ekin = 0.04 * std::fabs(normal());
0837 vecMass = targetParticle.GetMass() / GeV;
0838 totalEnergy = ekin + vecMass;
0839 targetParticle.SetTotalEnergy(totalEnergy * GeV);
0840 pp = std::sqrt(std::abs(totalEnergy * totalEnergy - vecMass * vecMass)) * GeV;
0841 pp1 = pseudoParticle[6].GetMomentum().mag() / MeV;
0842 if (pp1 < 1.0e-6 * GeV) {
0843 rthnve = pi * G4UniformRand();
0844 phinve = twopi * G4UniformRand();
0845 G4double srth = std::sin(rthnve);
0846 targetParticle.SetMomentum(
0847 pp * srth * std::cos(phinve) * MeV, pp * srth * std::sin(phinve) * MeV, pp * std::cos(rthnve) * MeV);
0848 } else {
0849 targetParticle.SetMomentum(pseudoParticle[6].GetMomentum() * (pp / pp1));
0850 }
0851 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
0852 } else
0853 {
0854 const G4double cpar[] = {0.6, 0.6, 0.35, 0.15, 0.10};
0855 const G4double gpar[] = {2.6, 2.6, 1.80, 1.30, 1.20};
0856
0857 G4int tempCount = std::max(1, std::min(5, backwardNucleonCount)) - 1;
0858
0859
0860 G4double rmb0 = 0.0;
0861 if (targetParticle.GetSide() == -3)
0862 rmb0 += targetParticle.GetMass() / GeV;
0863 for (i = 0; i < vecLen; ++i) {
0864 if (vec[i]->GetSide() == -3)
0865 rmb0 += vec[i]->GetMass() / GeV;
0866 }
0867 rmb = rmb0 + std::pow(-std::log(1.0 - G4UniformRand()), cpar[tempCount]) / gpar[tempCount];
0868 totalEnergy = pseudoParticle[6].GetTotalEnergy() / GeV;
0869 vecMass = std::min(rmb, totalEnergy);
0870 pseudoParticle[6].SetMass(vecMass * GeV);
0871 pp = std::sqrt(std::abs(totalEnergy * totalEnergy - vecMass * vecMass)) * GeV;
0872 pp1 = pseudoParticle[6].GetMomentum().mag() / MeV;
0873 if (pp1 < 1.0e-6 * GeV) {
0874 rthnve = pi * G4UniformRand();
0875 phinve = twopi * G4UniformRand();
0876 G4double srth = std::sin(rthnve);
0877 pseudoParticle[6].SetMomentum(
0878 -pp * srth * std::cos(phinve) * MeV, -pp * srth * std::sin(phinve) * MeV, -pp * std::cos(rthnve) * MeV);
0879 } else
0880 pseudoParticle[6].SetMomentum(pseudoParticle[6].GetMomentum() * (-pp / pp1));
0881
0882 G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV;
0883 tempV.Initialize(backwardNucleonCount);
0884 G4int tempLen = 0;
0885 if (targetParticle.GetSide() == -3)
0886 tempV.SetElement(tempLen++, &targetParticle);
0887 for (i = 0; i < vecLen; ++i) {
0888 if (vec[i]->GetSide() == -3)
0889 tempV.SetElement(tempLen++, vec[i]);
0890 }
0891 if (tempLen != backwardNucleonCount) {
0892 G4cerr << "tempLen is not the same as backwardNucleonCount" << G4endl;
0893 G4cerr << "tempLen = " << tempLen;
0894 G4cerr << ", backwardNucleonCount = " << backwardNucleonCount << G4endl;
0895 G4cerr << "targetParticle side = " << targetParticle.GetSide() << G4endl;
0896 G4cerr << "currentParticle side = " << currentParticle.GetSide() << G4endl;
0897 for (i = 0; i < vecLen; ++i)
0898 G4cerr << "particle #" << i << " side = " << vec[i]->GetSide() << G4endl;
0899 exit(EXIT_FAILURE);
0900 }
0901 constantCrossSection = true;
0902
0903 if (tempLen >= 2) {
0904 GenerateNBodyEvent(pseudoParticle[6].GetMass(), constantCrossSection, tempV, tempLen);
0905
0906 if (targetParticle.GetSide() == -3) {
0907 targetParticle.Lorentz(targetParticle, pseudoParticle[6]);
0908
0909 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
0910 }
0911 for (i = 0; i < vecLen; ++i) {
0912 if (vec[i]->GetSide() == -3) {
0913 vec[i]->Lorentz(*vec[i], pseudoParticle[6]);
0914 pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
0915 }
0916 }
0917
0918 }
0919 }
0920
0921
0922
0923 if (vecLen == 0)
0924 return false;
0925
0926
0927 G4int numberofFinalStateNucleons = 0;
0928 if (currentParticle.GetDefinition() == aProton || currentParticle.GetDefinition() == aNeutron ||
0929 currentParticle.GetDefinition() == G4SigmaMinus::SigmaMinus() ||
0930 currentParticle.GetDefinition() == G4SigmaPlus::SigmaPlus() ||
0931 currentParticle.GetDefinition() == G4SigmaZero::SigmaZero() ||
0932 currentParticle.GetDefinition() == G4XiZero::XiZero() ||
0933 currentParticle.GetDefinition() == G4XiMinus::XiMinus() ||
0934 currentParticle.GetDefinition() == G4OmegaMinus::OmegaMinus() ||
0935 currentParticle.GetDefinition() == G4Lambda::Lambda())
0936 ++numberofFinalStateNucleons;
0937 currentParticle.Lorentz(currentParticle, pseudoParticle[1]);
0938
0939 if (targetParticle.GetDefinition() == aProton || targetParticle.GetDefinition() == aNeutron ||
0940 targetParticle.GetDefinition() == G4Lambda::Lambda() || targetParticle.GetDefinition() == G4XiZero::XiZero() ||
0941 targetParticle.GetDefinition() == G4XiMinus::XiMinus() ||
0942 targetParticle.GetDefinition() == G4OmegaMinus::OmegaMinus() ||
0943 targetParticle.GetDefinition() == G4SigmaZero::SigmaZero() ||
0944 targetParticle.GetDefinition() == G4SigmaPlus::SigmaPlus() ||
0945 targetParticle.GetDefinition() == G4SigmaMinus::SigmaMinus())
0946 ++numberofFinalStateNucleons;
0947 if (targetParticle.GetDefinition() == G4AntiProton::AntiProton())
0948 --numberofFinalStateNucleons;
0949 if (targetParticle.GetDefinition() == G4AntiNeutron::AntiNeutron())
0950 --numberofFinalStateNucleons;
0951 if (targetParticle.GetDefinition() == G4AntiSigmaMinus::AntiSigmaMinus())
0952 --numberofFinalStateNucleons;
0953 if (targetParticle.GetDefinition() == G4AntiSigmaPlus::AntiSigmaPlus())
0954 --numberofFinalStateNucleons;
0955 if (targetParticle.GetDefinition() == G4AntiSigmaZero::AntiSigmaZero())
0956 --numberofFinalStateNucleons;
0957 if (targetParticle.GetDefinition() == G4AntiXiZero::AntiXiZero())
0958 --numberofFinalStateNucleons;
0959 if (targetParticle.GetDefinition() == G4AntiXiMinus::AntiXiMinus())
0960 --numberofFinalStateNucleons;
0961 if (targetParticle.GetDefinition() == G4AntiOmegaMinus::AntiOmegaMinus())
0962 --numberofFinalStateNucleons;
0963 if (targetParticle.GetDefinition() == G4AntiLambda::AntiLambda())
0964 --numberofFinalStateNucleons;
0965 targetParticle.Lorentz(targetParticle, pseudoParticle[1]);
0966
0967 for (i = 0; i < vecLen; ++i) {
0968 if (vec[i]->GetDefinition() == aProton || vec[i]->GetDefinition() == aNeutron ||
0969 vec[i]->GetDefinition() == G4Lambda::Lambda() || vec[i]->GetDefinition() == G4XiZero::XiZero() ||
0970 vec[i]->GetDefinition() == G4XiMinus::XiMinus() || vec[i]->GetDefinition() == G4OmegaMinus::OmegaMinus() ||
0971 vec[i]->GetDefinition() == G4SigmaPlus::SigmaPlus() || vec[i]->GetDefinition() == G4SigmaZero::SigmaZero() ||
0972 vec[i]->GetDefinition() == G4SigmaMinus::SigmaMinus())
0973 ++numberofFinalStateNucleons;
0974 if (vec[i]->GetDefinition() == G4AntiProton::AntiProton())
0975 --numberofFinalStateNucleons;
0976 if (vec[i]->GetDefinition() == G4AntiNeutron::AntiNeutron())
0977 --numberofFinalStateNucleons;
0978 if (vec[i]->GetDefinition() == G4AntiSigmaMinus::AntiSigmaMinus())
0979 --numberofFinalStateNucleons;
0980 if (vec[i]->GetDefinition() == G4AntiSigmaPlus::AntiSigmaPlus())
0981 --numberofFinalStateNucleons;
0982 if (vec[i]->GetDefinition() == G4AntiSigmaZero::AntiSigmaZero())
0983 --numberofFinalStateNucleons;
0984 if (vec[i]->GetDefinition() == G4AntiLambda::AntiLambda())
0985 --numberofFinalStateNucleons;
0986 if (vec[i]->GetDefinition() == G4AntiXiZero::AntiXiZero())
0987 --numberofFinalStateNucleons;
0988 if (vec[i]->GetDefinition() == G4AntiXiMinus::AntiXiMinus())
0989 --numberofFinalStateNucleons;
0990 if (vec[i]->GetDefinition() == G4AntiOmegaMinus::AntiOmegaMinus())
0991 --numberofFinalStateNucleons;
0992 vec[i]->Lorentz(*vec[i], pseudoParticle[1]);
0993 }
0994
0995 if (veryForward)
0996 numberofFinalStateNucleons++;
0997 numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
0998
0999
1000
1001
1002
1003
1004
1005
1006
1007
1008 G4bool leadingStrangeParticleHasChanged = true;
1009 if (leadFlag) {
1010 if (currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition())
1011 leadingStrangeParticleHasChanged = false;
1012 if (leadingStrangeParticleHasChanged && (targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition()))
1013 leadingStrangeParticleHasChanged = false;
1014 if (leadingStrangeParticleHasChanged) {
1015 for (i = 0; i < vecLen; i++) {
1016 if (vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition()) {
1017 leadingStrangeParticleHasChanged = false;
1018 break;
1019 }
1020 }
1021 }
1022 if (leadingStrangeParticleHasChanged) {
1023 G4bool leadTest =
1024 (leadingStrangeParticle.GetDefinition() == aKaonMinus ||
1025 leadingStrangeParticle.GetDefinition() == aKaonZeroL ||
1026 leadingStrangeParticle.GetDefinition() == aKaonZeroS ||
1027 leadingStrangeParticle.GetDefinition() == aKaonPlus || leadingStrangeParticle.GetDefinition() == aPiMinus ||
1028 leadingStrangeParticle.GetDefinition() == aPiZero || leadingStrangeParticle.GetDefinition() == aPiPlus);
1029 G4bool targetTest = false;
1030
1031
1032
1033 if ((leadTest && targetTest) || !(leadTest || targetTest))
1034 {
1035 targetParticle.SetDefinitionAndUpdateE(leadingStrangeParticle.GetDefinition());
1036 targetHasChanged = true;
1037
1038 } else {
1039 currentParticle.SetDefinitionAndUpdateE(leadingStrangeParticle.GetDefinition());
1040 incidentHasChanged = false;
1041
1042 }
1043 }
1044 }
1045
1046 pseudoParticle[3].SetMomentum(0.0, 0.0, pOriginal * GeV);
1047 pseudoParticle[3].SetMass(mOriginal * GeV);
1048 pseudoParticle[3].SetTotalEnergy(std::sqrt(pOriginal * pOriginal + mOriginal * mOriginal) * GeV);
1049
1050 const G4ParticleDefinition *aOrgDef = modifiedOriginal.GetDefinition();
1051 G4int diff = 0;
1052 if (aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron())
1053 diff = 1;
1054 if (numberofFinalStateNucleons == 1)
1055 diff = 0;
1056 pseudoParticle[4].SetMomentum(0.0, 0.0, 0.0);
1057 pseudoParticle[4].SetMass(protonMass * (numberofFinalStateNucleons - diff) * MeV);
1058 pseudoParticle[4].SetTotalEnergy(protonMass * (numberofFinalStateNucleons - diff) * MeV);
1059
1060 G4double theoreticalKinetic = pseudoParticle[3].GetTotalEnergy() / MeV + pseudoParticle[4].GetTotalEnergy() / MeV -
1061 currentParticle.GetMass() / MeV - targetParticle.GetMass() / MeV;
1062
1063 G4double simulatedKinetic = currentParticle.GetKineticEnergy() / MeV + targetParticle.GetKineticEnergy() / MeV;
1064
1065 pseudoParticle[5] = pseudoParticle[3] + pseudoParticle[4];
1066 pseudoParticle[3].Lorentz(pseudoParticle[3], pseudoParticle[5]);
1067 pseudoParticle[4].Lorentz(pseudoParticle[4], pseudoParticle[5]);
1068
1069 pseudoParticle[7].SetZero();
1070 pseudoParticle[7] = pseudoParticle[7] + currentParticle;
1071 pseudoParticle[7] = pseudoParticle[7] + targetParticle;
1072
1073 for (i = 0; i < vecLen; ++i) {
1074 pseudoParticle[7] = pseudoParticle[7] + *vec[i];
1075 simulatedKinetic += vec[i]->GetKineticEnergy() / MeV;
1076 theoreticalKinetic -= vec[i]->GetMass() / MeV;
1077 }
1078 if (vecLen <= 16 && vecLen > 0) {
1079
1080
1081
1082 G4ReactionProduct tempR[130];
1083 tempR[0] = currentParticle;
1084 tempR[1] = targetParticle;
1085 for (i = 0; i < vecLen; ++i)
1086 tempR[i + 2] = *vec[i];
1087 G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV;
1088 tempV.Initialize(vecLen + 2);
1089 G4int tempLen = 0;
1090 for (i = 0; i < vecLen + 2; ++i)
1091 tempV.SetElement(tempLen++, &tempR[i]);
1092 constantCrossSection = true;
1093
1094 wgt = GenerateNBodyEvent(pseudoParticle[3].GetTotalEnergy() / MeV + pseudoParticle[4].GetTotalEnergy() / MeV,
1095 constantCrossSection,
1096 tempV,
1097 tempLen);
1098 if (wgt > -.5) {
1099 theoreticalKinetic = 0.0;
1100 for (i = 0; i < tempLen; ++i) {
1101 pseudoParticle[6].Lorentz(*tempV[i], pseudoParticle[4]);
1102 theoreticalKinetic += pseudoParticle[6].GetKineticEnergy() / MeV;
1103 }
1104 }
1105
1106 }
1107
1108
1109
1110 if (simulatedKinetic != 0.0) {
1111 wgt = (theoreticalKinetic) / simulatedKinetic;
1112 theoreticalKinetic = currentParticle.GetKineticEnergy() / MeV * wgt;
1113 simulatedKinetic = theoreticalKinetic;
1114 currentParticle.SetKineticEnergy(theoreticalKinetic * MeV);
1115 pp = currentParticle.GetTotalMomentum() / MeV;
1116 pp1 = currentParticle.GetMomentum().mag() / MeV;
1117 if (pp1 < 1.0e-6 * GeV) {
1118 rthnve = pi * G4UniformRand();
1119 phinve = twopi * G4UniformRand();
1120 currentParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1121 pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1122 pp * std::cos(rthnve) * MeV);
1123 } else {
1124 currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
1125 }
1126 theoreticalKinetic = targetParticle.GetKineticEnergy() / MeV * wgt;
1127 targetParticle.SetKineticEnergy(theoreticalKinetic * MeV);
1128 simulatedKinetic += theoreticalKinetic;
1129 pp = targetParticle.GetTotalMomentum() / MeV;
1130 pp1 = targetParticle.GetMomentum().mag() / MeV;
1131
1132 if (pp1 < 1.0e-6 * GeV) {
1133 rthnve = pi * G4UniformRand();
1134 phinve = twopi * G4UniformRand();
1135 targetParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1136 pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1137 pp * std::cos(rthnve) * MeV);
1138 } else {
1139 targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
1140 }
1141 for (i = 0; i < vecLen; ++i) {
1142 theoreticalKinetic = vec[i]->GetKineticEnergy() / MeV * wgt;
1143 simulatedKinetic += theoreticalKinetic;
1144 vec[i]->SetKineticEnergy(theoreticalKinetic * MeV);
1145 pp = vec[i]->GetTotalMomentum() / MeV;
1146 pp1 = vec[i]->GetMomentum().mag() / MeV;
1147 if (pp1 < 1.0e-6 * GeV) {
1148 rthnve = pi * G4UniformRand();
1149 phinve = twopi * G4UniformRand();
1150 vec[i]->SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1151 pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1152 pp * std::cos(rthnve) * MeV);
1153 } else
1154 vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp / pp1));
1155 }
1156 }
1157
1158 Rotate(numberofFinalStateNucleons,
1159 pseudoParticle[3].GetMomentum(),
1160 modifiedOriginal,
1161 originalIncident,
1162 targetNucleus,
1163 currentParticle,
1164 targetParticle,
1165 vec,
1166 vecLen);
1167
1168
1169
1170
1171
1172
1173 if (atomicWeight >= 1.5) {
1174
1175
1176
1177
1178
1179 G4double epnb, edta;
1180 G4int npnb = 0;
1181 G4int ndta = 0;
1182
1183 epnb = targetNucleus.GetPNBlackTrackEnergy();
1184 edta = targetNucleus.GetDTABlackTrackEnergy();
1185 const G4double pnCutOff = 0.001;
1186 const G4double dtaCutOff = 0.001;
1187 const G4double kineticMinimum = 1.e-6;
1188 const G4double kineticFactor = -0.010;
1189 G4double sprob = 0.0;
1190 const G4double ekIncident = originalIncident->GetKineticEnergy() / GeV;
1191 if (ekIncident >= 5.0)
1192 sprob = std::min(1.0, 0.6 * std::log(ekIncident - 4.0));
1193 if (epnb >= pnCutOff) {
1194 npnb = Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * epnb / (epnb + edta));
1195 if (numberofFinalStateNucleons + npnb > atomicWeight)
1196 npnb = G4int(atomicWeight + 0.00001 - numberofFinalStateNucleons);
1197 npnb = std::min(npnb, 127 - vecLen);
1198 }
1199 if (edta >= dtaCutOff) {
1200 ndta = Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * edta / (epnb + edta));
1201 ndta = std::min(ndta, 127 - vecLen);
1202 }
1203 G4double spall = numberofFinalStateNucleons;
1204
1205
1206 AddBlackTrackParticles(epnb,
1207 npnb,
1208 edta,
1209 ndta,
1210 sprob,
1211 kineticMinimum,
1212 kineticFactor,
1213 modifiedOriginal,
1214 spall,
1215 targetNucleus,
1216 vec,
1217 vecLen);
1218
1219
1220 }
1221
1222
1223
1224 if ((atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2))
1225 currentParticle.SetTOF(1.0 - 500.0 * std::exp(-ekOriginal / 0.04) * std::log(G4UniformRand()));
1226 else
1227 currentParticle.SetTOF(1.0);
1228 return true;
1229
1230 }
1231
1232 void FullModelReactionDynamics::SuppressChargedPions(G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
1233 G4int &vecLen,
1234 const G4ReactionProduct &modifiedOriginal,
1235 G4ReactionProduct ¤tParticle,
1236 G4ReactionProduct &targetParticle,
1237 const G4Nucleus &targetNucleus,
1238 G4bool &incidentHasChanged,
1239 G4bool &targetHasChanged) {
1240
1241
1242
1243
1244 const G4double atomicWeight = targetNucleus.GetN_asInt();
1245 const G4double atomicNumber = targetNucleus.GetZ_asInt();
1246 const G4double pOriginal = modifiedOriginal.GetTotalMomentum() / GeV;
1247
1248
1249 G4ParticleDefinition *aProton = G4Proton::Proton();
1250 G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
1251 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
1252 G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
1253 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
1254 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
1255
1256 const G4bool antiTest =
1257 modifiedOriginal.GetDefinition() != anAntiProton && modifiedOriginal.GetDefinition() != anAntiNeutron;
1258 if (antiTest &&
1259 (
1260
1261 currentParticle.GetDefinition() == aPiPlus || currentParticle.GetDefinition() == aPiMinus) &&
1262 (G4UniformRand() <= (10.0 - pOriginal) / 6.0) && (G4UniformRand() <= atomicWeight / 300.0)) {
1263 if (G4UniformRand() > atomicNumber / atomicWeight)
1264 currentParticle.SetDefinitionAndUpdateE(aNeutron);
1265 else
1266 currentParticle.SetDefinitionAndUpdateE(aProton);
1267 incidentHasChanged = true;
1268 }
1269 for (G4int i = 0; i < vecLen; ++i) {
1270 if (antiTest &&
1271 (
1272
1273 vec[i]->GetDefinition() == aPiPlus || vec[i]->GetDefinition() == aPiMinus) &&
1274 (G4UniformRand() <= (10.0 - pOriginal) / 6.0) && (G4UniformRand() <= atomicWeight / 300.0)) {
1275 if (G4UniformRand() > atomicNumber / atomicWeight)
1276 vec[i]->SetDefinitionAndUpdateE(aNeutron);
1277 else
1278 vec[i]->SetDefinitionAndUpdateE(aProton);
1279
1280 }
1281 }
1282
1283 }
1284
1285 G4bool FullModelReactionDynamics::TwoCluster(
1286 G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
1287 G4int &vecLen,
1288 G4ReactionProduct &modifiedOriginal,
1289 const G4HadProjectile *originalIncident,
1290 G4ReactionProduct ¤tParticle,
1291 G4ReactionProduct &targetParticle,
1292 const G4Nucleus &targetNucleus,
1293 G4bool &incidentHasChanged,
1294 G4bool &targetHasChanged,
1295 G4bool leadFlag,
1296 G4ReactionProduct &leadingStrangeParticle) {
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310 G4int i;
1311 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
1312 G4ParticleDefinition *aProton = G4Proton::Proton();
1313 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
1314 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
1315 G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
1316 G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
1317 const G4double protonMass = aProton->GetPDGMass() / MeV;
1318 const G4double ekOriginal = modifiedOriginal.GetKineticEnergy() / GeV;
1319 const G4double etOriginal = modifiedOriginal.GetTotalEnergy() / GeV;
1320 const G4double mOriginal = modifiedOriginal.GetMass() / GeV;
1321 const G4double pOriginal = modifiedOriginal.GetMomentum().mag() / GeV;
1322 G4double targetMass = targetParticle.GetDefinition()->GetPDGMass() / GeV;
1323 G4double centerofmassEnergy =
1324 std::sqrt(mOriginal * mOriginal + targetMass * targetMass + 2.0 * targetMass * etOriginal);
1325 G4double currentMass = currentParticle.GetMass() / GeV;
1326 targetMass = targetParticle.GetMass() / GeV;
1327
1328 if (currentMass == 0.0 && targetMass == 0.0) {
1329 G4double ek = currentParticle.GetKineticEnergy();
1330 G4ThreeVector m = currentParticle.GetMomentum();
1331 currentParticle = *vec[0];
1332 targetParticle = *vec[1];
1333 for (i = 0; i < (vecLen - 2); ++i)
1334 *vec[i] = *vec[i + 2];
1335 if (vecLen < 2) {
1336 for (G4int i = 0; i < vecLen; i++)
1337 delete vec[i];
1338 vecLen = 0;
1339 G4ExceptionDescription ed;
1340 ed << "Negative number of particles";
1341 G4Exception("FullModelReactionDynamics::TwoCluster", "had064", FatalException, ed);
1342 }
1343 delete vec[vecLen - 1];
1344 delete vec[vecLen - 2];
1345 vecLen -= 2;
1346 incidentHasChanged = true;
1347 targetHasChanged = true;
1348 currentParticle.SetKineticEnergy(ek);
1349 currentParticle.SetMomentum(m);
1350 }
1351 const G4double atomicWeight = targetNucleus.GetN_asInt();
1352 const G4double atomicNumber = targetNucleus.GetZ_asInt();
1353
1354
1355
1356
1357
1358 G4int forwardCount = 1;
1359 currentParticle.SetSide(1);
1360 G4double forwardMass = currentParticle.GetMass() / GeV;
1361
1362
1363
1364 G4int backwardCount = 1;
1365 G4int backwardNucleonCount = 1;
1366 targetParticle.SetSide(-1);
1367 G4double backwardMass = targetParticle.GetMass() / GeV;
1368
1369
1370 for (i = 0; i < vecLen; ++i) {
1371 if (vec[i]->GetSide() < 0)
1372 vec[i]->SetSide(-1);
1373
1374
1375 if (vec[i]->GetSide() == -1) {
1376 ++backwardCount;
1377 backwardMass += vec[i]->GetMass() / GeV;
1378 } else {
1379 ++forwardCount;
1380 forwardMass += vec[i]->GetMass() / GeV;
1381 }
1382 }
1383
1384
1385
1386 G4double term1 = std::log(centerofmassEnergy * centerofmassEnergy);
1387 if (term1 < 0)
1388 term1 = 0.0001;
1389 const G4double afc = 0.312 + 0.2 * std::log(term1);
1390 G4double xtarg;
1391 if (centerofmassEnergy < 2.0 + G4UniformRand())
1392 xtarg = afc * (std::pow(atomicWeight, 0.33) - 1.0) * (2 * backwardCount + vecLen + 2) / 2.0;
1393 else
1394 xtarg = afc * (std::pow(atomicWeight, 0.33) - 1.0) * (2 * backwardCount);
1395 if (xtarg <= 0.0)
1396 xtarg = 0.01;
1397 G4int nuclearExcitationCount = Poisson(xtarg);
1398 if (atomicWeight < 1.0001)
1399 nuclearExcitationCount = 0;
1400 G4int extraNucleonCount = 0;
1401
1402
1403 if (nuclearExcitationCount > 0) {
1404 G4int momentumBin = std::min(4, G4int(pOriginal / 3.0));
1405 const G4double nucsup[] = {1.0, 0.8, 0.6, 0.5, 0.4};
1406
1407
1408
1409
1410 for (i = 0; i < nuclearExcitationCount; ++i) {
1411 G4ReactionProduct *pVec = new G4ReactionProduct();
1412 if (G4UniformRand() < nucsup[momentumBin])
1413 {
1414 if (G4UniformRand() > 1.0 - atomicNumber / atomicWeight)
1415
1416 pVec->SetDefinition(aProton);
1417 else
1418 pVec->SetDefinition(aNeutron);
1419 ++backwardNucleonCount;
1420 ++extraNucleonCount;
1421
1422 } else {
1423 G4double ran = G4UniformRand();
1424 if (ran < 0.3181)
1425 pVec->SetDefinition(aPiPlus);
1426 else if (ran < 0.6819)
1427 pVec->SetDefinition(aPiZero);
1428 else
1429 pVec->SetDefinition(aPiMinus);
1430 }
1431 pVec->SetSide(-2);
1432
1433 pVec->SetNewlyAdded(true);
1434 vec.SetElement(vecLen++, pVec);
1435
1436 }
1437 }
1438
1439
1440
1441 G4double eAvailable = centerofmassEnergy - (forwardMass + backwardMass);
1442 G4bool secondaryDeleted;
1443 G4double pMass;
1444 while (eAvailable <= 0.0)
1445 {
1446 secondaryDeleted = false;
1447 for (i = (vecLen - 1); i >= 0; --i) {
1448 if (vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled()) {
1449 pMass = vec[i]->GetMass() / GeV;
1450 for (G4int j = i; j < (vecLen - 1); ++j)
1451 *vec[j] = *vec[j + 1];
1452 --forwardCount;
1453
1454 forwardMass -= pMass;
1455 secondaryDeleted = true;
1456 break;
1457 } else if (vec[i]->GetSide() == -1 && vec[i]->GetMayBeKilled()) {
1458 pMass = vec[i]->GetMass() / GeV;
1459 for (G4int j = i; j < (vecLen - 1); ++j)
1460 *vec[j] = *vec[j + 1];
1461 --backwardCount;
1462
1463 backwardMass -= pMass;
1464 secondaryDeleted = true;
1465 break;
1466 }
1467 }
1468 if (secondaryDeleted) {
1469 G4ReactionProduct *temp = vec[vecLen - 1];
1470 delete temp;
1471 --vecLen;
1472
1473 } else {
1474 if (vecLen == 0) {
1475 return false;
1476 }
1477 if (targetParticle.GetSide() == -1) {
1478 pMass = targetParticle.GetMass() / GeV;
1479 targetParticle = *vec[0];
1480 for (G4int j = 0; j < (vecLen - 1); ++j)
1481 *vec[j] = *vec[j + 1];
1482 --backwardCount;
1483
1484 backwardMass -= pMass;
1485 secondaryDeleted = true;
1486 } else if (targetParticle.GetSide() == 1) {
1487 pMass = targetParticle.GetMass() / GeV;
1488 targetParticle = *vec[0];
1489 for (G4int j = 0; j < (vecLen - 1); ++j)
1490 *vec[j] = *vec[j + 1];
1491 --forwardCount;
1492
1493 forwardMass -= pMass;
1494 secondaryDeleted = true;
1495 }
1496 if (secondaryDeleted) {
1497 G4ReactionProduct *temp = vec[vecLen - 1];
1498 delete temp;
1499 --vecLen;
1500
1501 } else {
1502 if (currentParticle.GetSide() == -1) {
1503 pMass = currentParticle.GetMass() / GeV;
1504 currentParticle = *vec[0];
1505 for (G4int j = 0; j < (vecLen - 1); ++j)
1506 *vec[j] = *vec[j + 1];
1507 --backwardCount;
1508
1509 backwardMass -= pMass;
1510 secondaryDeleted = true;
1511 } else if (currentParticle.GetSide() == 1) {
1512 pMass = currentParticle.GetMass() / GeV;
1513 currentParticle = *vec[0];
1514 for (G4int j = 0; j < (vecLen - 1); ++j)
1515 *vec[j] = *vec[j + 1];
1516 --forwardCount;
1517
1518 forwardMass -= pMass;
1519 secondaryDeleted = true;
1520 }
1521 if (secondaryDeleted) {
1522 G4ReactionProduct *temp = vec[vecLen - 1];
1523 delete temp;
1524 --vecLen;
1525
1526 } else
1527 break;
1528 }
1529 }
1530 eAvailable = centerofmassEnergy - (forwardMass + backwardMass);
1531 }
1532
1533
1534
1535
1536
1537
1538
1539 G4double rmc = 0.0, rmd = 0.0;
1540 const G4double cpar[] = {0.6, 0.6, 0.35, 0.15, 0.10};
1541 const G4double gpar[] = {2.6, 2.6, 1.8, 1.30, 1.20};
1542
1543 if (forwardCount == 0)
1544 return false;
1545
1546 if (forwardCount == 1)
1547 rmc = forwardMass;
1548 else {
1549 G4int ntc = std::max(1, std::min(5, forwardCount)) - 1;
1550 rmc = forwardMass + std::pow(-std::log(1.0 - G4UniformRand()), cpar[ntc]) / gpar[ntc];
1551 }
1552 if (backwardCount == 1)
1553 rmd = backwardMass;
1554 else {
1555 G4int ntc = std::max(1, std::min(5, backwardCount)) - 1;
1556 rmd = backwardMass + std::pow(-std::log(1.0 - G4UniformRand()), cpar[ntc]) / gpar[ntc];
1557 }
1558 while (rmc + rmd > centerofmassEnergy) {
1559 if ((rmc <= forwardMass) && (rmd <= backwardMass)) {
1560 G4double temp = 0.999 * centerofmassEnergy / (rmc + rmd);
1561 rmc *= temp;
1562 rmd *= temp;
1563 } else {
1564 rmc = 0.1 * forwardMass + 0.9 * rmc;
1565 rmd = 0.1 * backwardMass + 0.9 * rmd;
1566 }
1567 }
1568
1569
1570
1571 G4ReactionProduct pseudoParticle[8];
1572 for (i = 0; i < 8; ++i)
1573 pseudoParticle[i].SetZero();
1574
1575 pseudoParticle[1].SetMass(mOriginal * GeV);
1576 pseudoParticle[1].SetTotalEnergy(etOriginal * GeV);
1577 pseudoParticle[1].SetMomentum(0.0, 0.0, pOriginal * GeV);
1578
1579 pseudoParticle[2].SetMass(protonMass * MeV);
1580 pseudoParticle[2].SetTotalEnergy(protonMass * MeV);
1581 pseudoParticle[2].SetMomentum(0.0, 0.0, 0.0);
1582
1583
1584
1585 pseudoParticle[0] = pseudoParticle[1] + pseudoParticle[2];
1586 pseudoParticle[1].Lorentz(pseudoParticle[1], pseudoParticle[0]);
1587 pseudoParticle[2].Lorentz(pseudoParticle[2], pseudoParticle[0]);
1588
1589 const G4double pfMin = 0.0001;
1590 G4double pf = (centerofmassEnergy * centerofmassEnergy + rmd * rmd - rmc * rmc);
1591 pf *= pf;
1592 pf -= 4 * centerofmassEnergy * centerofmassEnergy * rmd * rmd;
1593 pf = std::sqrt(std::max(pf, pfMin)) / (2.0 * centerofmassEnergy);
1594
1595
1596
1597 pseudoParticle[3].SetMass(rmc * GeV);
1598 pseudoParticle[3].SetTotalEnergy(std::sqrt(pf * pf + rmc * rmc) * GeV);
1599
1600 pseudoParticle[4].SetMass(rmd * GeV);
1601 pseudoParticle[4].SetTotalEnergy(std::sqrt(pf * pf + rmd * rmd) * GeV);
1602
1603
1604
1605 const G4double bMin = 0.01;
1606 const G4double b1 = 4.0;
1607 const G4double b2 = 1.6;
1608 G4double t = std::log(1.0 - G4UniformRand()) / std::max(bMin, b1 + b2 * std::log(pOriginal));
1609 G4double t1 = pseudoParticle[1].GetTotalEnergy() / GeV - pseudoParticle[3].GetTotalEnergy() / GeV;
1610 G4double pin = pseudoParticle[1].GetMomentum().mag() / GeV;
1611 G4double tacmin = t1 * t1 - (pin - pf) * (pin - pf);
1612
1613
1614
1615 const G4double smallValue = 1.0e-10;
1616 G4double dumnve = 4.0 * pin * pf;
1617 if (dumnve == 0.0)
1618 dumnve = smallValue;
1619 G4double ctet = std::max(-1.0, std::min(1.0, 1.0 + 2.0 * (t - tacmin) / dumnve));
1620 dumnve = std::max(0.0, 1.0 - ctet * ctet);
1621 G4double stet = std::sqrt(dumnve);
1622 G4double phi = G4UniformRand() * twopi;
1623
1624
1625
1626 pseudoParticle[3].SetMomentum(pf * stet * std::sin(phi) * GeV, pf * stet * std::cos(phi) * GeV, pf * ctet * GeV);
1627 pseudoParticle[4].SetMomentum(pseudoParticle[3].GetMomentum() * (-1.0));
1628
1629
1630
1631 G4double pp, pp1, rthnve, phinve;
1632 if (nuclearExcitationCount > 0) {
1633 const G4double ga = 1.2;
1634 G4double ekit1 = 0.04;
1635 G4double ekit2 = 0.6;
1636 if (ekOriginal <= 5.0) {
1637 ekit1 *= ekOriginal * ekOriginal / 25.0;
1638 ekit2 *= ekOriginal * ekOriginal / 25.0;
1639 }
1640 const G4double a = (1.0 - ga) / (std::pow(ekit2, (1.0 - ga)) - std::pow(ekit1, (1.0 - ga)));
1641 for (i = 0; i < vecLen; ++i) {
1642 if (vec[i]->GetSide() == -2) {
1643 G4double kineticE =
1644 std::pow((G4UniformRand() * (1.0 - ga) / a + std::pow(ekit1, (1.0 - ga))), (1.0 / (1.0 - ga)));
1645 vec[i]->SetKineticEnergy(kineticE * GeV);
1646 G4double vMass = vec[i]->GetMass() / MeV;
1647 G4double totalE = kineticE + vMass;
1648 pp = std::sqrt(std::abs(totalE * totalE - vMass * vMass));
1649 G4double cost = std::min(1.0, std::max(-1.0, std::log(2.23 * G4UniformRand() + 0.383) / 0.96));
1650 G4double sint = std::sqrt(std::max(0.0, (1.0 - cost * cost)));
1651 phi = twopi * G4UniformRand();
1652 vec[i]->SetMomentum(pp * sint * std::sin(phi) * MeV, pp * sint * std::cos(phi) * MeV, pp * cost * MeV);
1653 vec[i]->Lorentz(*vec[i], pseudoParticle[0]);
1654 }
1655 }
1656
1657 }
1658
1659
1660
1661 currentParticle.SetMomentum(pseudoParticle[3].GetMomentum());
1662 currentParticle.SetTotalEnergy(pseudoParticle[3].GetTotalEnergy());
1663
1664 targetParticle.SetMomentum(pseudoParticle[4].GetMomentum());
1665 targetParticle.SetTotalEnergy(pseudoParticle[4].GetTotalEnergy());
1666
1667 pseudoParticle[5].SetMomentum(pseudoParticle[3].GetMomentum() * (-1.0));
1668 pseudoParticle[5].SetMass(pseudoParticle[3].GetMass());
1669 pseudoParticle[5].SetTotalEnergy(pseudoParticle[3].GetTotalEnergy());
1670
1671 pseudoParticle[6].SetMomentum(pseudoParticle[4].GetMomentum() * (-1.0));
1672 pseudoParticle[6].SetMass(pseudoParticle[4].GetMass());
1673 pseudoParticle[6].SetTotalEnergy(pseudoParticle[4].GetTotalEnergy());
1674
1675 G4double wgt;
1676
1677 if (forwardCount > 1)
1678 {
1679 G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV;
1680 tempV.Initialize(forwardCount);
1681 G4bool constantCrossSection = true;
1682 G4int tempLen = 0;
1683 if (currentParticle.GetSide() == 1)
1684 tempV.SetElement(tempLen++, ¤tParticle);
1685 if (targetParticle.GetSide() == 1)
1686 tempV.SetElement(tempLen++, &targetParticle);
1687 for (i = 0; i < vecLen; ++i) {
1688 if (vec[i]->GetSide() == 1) {
1689 if (tempLen < 18)
1690 tempV.SetElement(tempLen++, vec[i]);
1691 else {
1692 vec[i]->SetSide(-1);
1693 continue;
1694 }
1695 }
1696 }
1697 if (tempLen >= 2) {
1698 GenerateNBodyEvent(pseudoParticle[3].GetMass() / MeV, constantCrossSection, tempV, tempLen);
1699 if (currentParticle.GetSide() == 1)
1700 currentParticle.Lorentz(currentParticle, pseudoParticle[5]);
1701 if (targetParticle.GetSide() == 1)
1702 targetParticle.Lorentz(targetParticle, pseudoParticle[5]);
1703 for (i = 0; i < vecLen; ++i) {
1704 if (vec[i]->GetSide() == 1)
1705 vec[i]->Lorentz(*vec[i], pseudoParticle[5]);
1706 }
1707 }
1708 }
1709
1710 if (backwardCount > 1)
1711 {
1712 G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV;
1713 tempV.Initialize(backwardCount);
1714 G4bool constantCrossSection = true;
1715 G4int tempLen = 0;
1716 if (currentParticle.GetSide() == -1)
1717 tempV.SetElement(tempLen++, ¤tParticle);
1718 if (targetParticle.GetSide() == -1)
1719 tempV.SetElement(tempLen++, &targetParticle);
1720 for (i = 0; i < vecLen; ++i) {
1721 if (vec[i]->GetSide() == -1) {
1722 if (tempLen < 18)
1723 tempV.SetElement(tempLen++, vec[i]);
1724 else {
1725 vec[i]->SetSide(-2);
1726 vec[i]->SetKineticEnergy(0.0);
1727 vec[i]->SetMomentum(0.0, 0.0, 0.0);
1728 continue;
1729 }
1730 }
1731 }
1732 if (tempLen >= 2) {
1733 GenerateNBodyEvent(pseudoParticle[4].GetMass() / MeV, constantCrossSection, tempV, tempLen);
1734 if (currentParticle.GetSide() == -1)
1735 currentParticle.Lorentz(currentParticle, pseudoParticle[6]);
1736 if (targetParticle.GetSide() == -1)
1737 targetParticle.Lorentz(targetParticle, pseudoParticle[6]);
1738 for (i = 0; i < vecLen; ++i) {
1739 if (vec[i]->GetSide() == -1)
1740 vec[i]->Lorentz(*vec[i], pseudoParticle[6]);
1741 }
1742 }
1743 }
1744
1745
1746
1747
1748 G4int numberofFinalStateNucleons = 0;
1749 if (currentParticle.GetDefinition() == aProton || currentParticle.GetDefinition() == aNeutron ||
1750 currentParticle.GetDefinition() == aSigmaMinus || currentParticle.GetDefinition() == G4SigmaPlus::SigmaPlus() ||
1751 currentParticle.GetDefinition() == G4SigmaZero::SigmaZero() ||
1752 currentParticle.GetDefinition() == G4XiZero::XiZero() ||
1753 currentParticle.GetDefinition() == G4XiMinus::XiMinus() ||
1754 currentParticle.GetDefinition() == G4OmegaMinus::OmegaMinus() ||
1755 currentParticle.GetDefinition() == G4Lambda::Lambda())
1756 ++numberofFinalStateNucleons;
1757 currentParticle.Lorentz(currentParticle, pseudoParticle[2]);
1758 if (targetParticle.GetDefinition() == aProton || targetParticle.GetDefinition() == aNeutron ||
1759 targetParticle.GetDefinition() == G4Lambda::Lambda() || targetParticle.GetDefinition() == G4XiZero::XiZero() ||
1760 targetParticle.GetDefinition() == G4XiMinus::XiMinus() ||
1761 targetParticle.GetDefinition() == G4OmegaMinus::OmegaMinus() ||
1762 targetParticle.GetDefinition() == G4SigmaZero::SigmaZero() ||
1763 targetParticle.GetDefinition() == G4SigmaPlus::SigmaPlus() || targetParticle.GetDefinition() == aSigmaMinus)
1764 ++numberofFinalStateNucleons;
1765 if (targetParticle.GetDefinition() == G4AntiProton::AntiProton())
1766 --numberofFinalStateNucleons;
1767 if (targetParticle.GetDefinition() == G4AntiNeutron::AntiNeutron())
1768 --numberofFinalStateNucleons;
1769 if (targetParticle.GetDefinition() == G4AntiSigmaMinus::AntiSigmaMinus())
1770 --numberofFinalStateNucleons;
1771 if (targetParticle.GetDefinition() == G4AntiSigmaPlus::AntiSigmaPlus())
1772 --numberofFinalStateNucleons;
1773 if (targetParticle.GetDefinition() == G4AntiSigmaZero::AntiSigmaZero())
1774 --numberofFinalStateNucleons;
1775 if (targetParticle.GetDefinition() == G4AntiXiZero::AntiXiZero())
1776 --numberofFinalStateNucleons;
1777 if (targetParticle.GetDefinition() == G4AntiXiMinus::AntiXiMinus())
1778 --numberofFinalStateNucleons;
1779 if (targetParticle.GetDefinition() == G4AntiOmegaMinus::AntiOmegaMinus())
1780 --numberofFinalStateNucleons;
1781 if (targetParticle.GetDefinition() == G4AntiLambda::AntiLambda())
1782 --numberofFinalStateNucleons;
1783 targetParticle.Lorentz(targetParticle, pseudoParticle[2]);
1784 for (i = 0; i < vecLen; ++i) {
1785 if (vec[i]->GetDefinition() == aProton || vec[i]->GetDefinition() == aNeutron ||
1786 vec[i]->GetDefinition() == G4Lambda::Lambda() || vec[i]->GetDefinition() == G4XiZero::XiZero() ||
1787 vec[i]->GetDefinition() == G4XiMinus::XiMinus() || vec[i]->GetDefinition() == G4OmegaMinus::OmegaMinus() ||
1788 vec[i]->GetDefinition() == G4SigmaPlus::SigmaPlus() || vec[i]->GetDefinition() == G4SigmaZero::SigmaZero() ||
1789 vec[i]->GetDefinition() == aSigmaMinus)
1790 ++numberofFinalStateNucleons;
1791 if (vec[i]->GetDefinition() == G4AntiProton::AntiProton())
1792 --numberofFinalStateNucleons;
1793 if (vec[i]->GetDefinition() == G4AntiNeutron::AntiNeutron())
1794 --numberofFinalStateNucleons;
1795 if (vec[i]->GetDefinition() == G4AntiSigmaMinus::AntiSigmaMinus())
1796 --numberofFinalStateNucleons;
1797 if (vec[i]->GetDefinition() == G4AntiSigmaPlus::AntiSigmaPlus())
1798 --numberofFinalStateNucleons;
1799 if (vec[i]->GetDefinition() == G4AntiSigmaZero::AntiSigmaZero())
1800 --numberofFinalStateNucleons;
1801 if (vec[i]->GetDefinition() == G4AntiLambda::AntiLambda())
1802 --numberofFinalStateNucleons;
1803 if (vec[i]->GetDefinition() == G4AntiXiZero::AntiXiZero())
1804 --numberofFinalStateNucleons;
1805 if (vec[i]->GetDefinition() == G4AntiXiMinus::AntiXiMinus())
1806 --numberofFinalStateNucleons;
1807 if (vec[i]->GetDefinition() == G4AntiOmegaMinus::AntiOmegaMinus())
1808 --numberofFinalStateNucleons;
1809 vec[i]->Lorentz(*vec[i], pseudoParticle[2]);
1810 }
1811
1812 numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
1813
1814
1815
1816 G4bool dum = true;
1817 if (leadFlag) {
1818
1819
1820
1821
1822
1823
1824
1825
1826
1827 if (currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition())
1828 dum = false;
1829 else if (targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition())
1830 dum = false;
1831 else {
1832 for (i = 0; i < vecLen; ++i) {
1833 if (vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition()) {
1834 dum = false;
1835 break;
1836 }
1837 }
1838 }
1839 if (dum) {
1840 G4double leadMass = leadingStrangeParticle.GetMass() / MeV;
1841 G4double ekin;
1842 if (((leadMass < protonMass) && (targetParticle.GetMass() / MeV < protonMass)) ||
1843 ((leadMass >= protonMass) && (targetParticle.GetMass() / MeV >= protonMass))) {
1844 ekin = targetParticle.GetKineticEnergy() / GeV;
1845 pp1 = targetParticle.GetMomentum().mag() / MeV;
1846 targetParticle.SetDefinition(leadingStrangeParticle.GetDefinition());
1847 targetParticle.SetKineticEnergy(ekin * GeV);
1848 pp = targetParticle.GetTotalMomentum() / MeV;
1849 if (pp1 < 1.0e-3) {
1850 rthnve = pi * G4UniformRand();
1851 phinve = twopi * G4UniformRand();
1852 targetParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1853 pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1854 pp * std::cos(rthnve) * MeV);
1855 } else
1856 targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
1857
1858 targetHasChanged = true;
1859 } else {
1860 ekin = currentParticle.GetKineticEnergy() / GeV;
1861 pp1 = currentParticle.GetMomentum().mag() / MeV;
1862 currentParticle.SetDefinition(leadingStrangeParticle.GetDefinition());
1863 currentParticle.SetKineticEnergy(ekin * GeV);
1864 pp = currentParticle.GetTotalMomentum() / MeV;
1865 if (pp1 < 1.0e-3) {
1866 rthnve = pi * G4UniformRand();
1867 phinve = twopi * G4UniformRand();
1868 currentParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1869 pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1870 pp * std::cos(rthnve) * MeV);
1871 } else
1872 currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
1873
1874 incidentHasChanged = true;
1875 }
1876 }
1877 }
1878
1879
1880
1881
1882 pseudoParticle[4].SetMass(mOriginal * GeV);
1883 pseudoParticle[4].SetTotalEnergy(etOriginal * GeV);
1884 pseudoParticle[4].SetMomentum(0.0, 0.0, pOriginal * GeV);
1885
1886 const G4ParticleDefinition *aOrgDef = modifiedOriginal.GetDefinition();
1887 G4int diff = 0;
1888 if (aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron())
1889 diff = 1;
1890 if (numberofFinalStateNucleons == 1)
1891 diff = 0;
1892 pseudoParticle[5].SetMomentum(0.0, 0.0, 0.0);
1893 pseudoParticle[5].SetMass(protonMass * (numberofFinalStateNucleons - diff) * MeV);
1894 pseudoParticle[5].SetTotalEnergy(protonMass * (numberofFinalStateNucleons - diff) * MeV);
1895
1896
1897 G4double theoreticalKinetic = pseudoParticle[4].GetTotalEnergy() / GeV + pseudoParticle[5].GetTotalEnergy() / GeV;
1898
1899 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
1900 pseudoParticle[4].Lorentz(pseudoParticle[4], pseudoParticle[6]);
1901 pseudoParticle[5].Lorentz(pseudoParticle[5], pseudoParticle[6]);
1902
1903 if (vecLen < 16) {
1904 G4ReactionProduct tempR[130];
1905
1906 tempR[0] = currentParticle;
1907 tempR[1] = targetParticle;
1908 for (i = 0; i < vecLen; ++i)
1909 tempR[i + 2] = *vec[i];
1910
1911 G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV;
1912 tempV.Initialize(vecLen + 2);
1913 G4bool constantCrossSection = true;
1914 G4int tempLen = 0;
1915 for (i = 0; i < vecLen + 2; ++i)
1916 tempV.SetElement(tempLen++, &tempR[i]);
1917
1918 if (tempLen >= 2) {
1919
1920 GenerateNBodyEvent(pseudoParticle[4].GetTotalEnergy() / MeV + pseudoParticle[5].GetTotalEnergy() / MeV,
1921 constantCrossSection,
1922 tempV,
1923 tempLen);
1924 theoreticalKinetic = 0.0;
1925 for (i = 0; i < vecLen + 2; ++i) {
1926 pseudoParticle[7].SetMomentum(tempV[i]->GetMomentum());
1927 pseudoParticle[7].SetMass(tempV[i]->GetMass());
1928 pseudoParticle[7].SetTotalEnergy(tempV[i]->GetTotalEnergy());
1929 pseudoParticle[7].Lorentz(pseudoParticle[7], pseudoParticle[5]);
1930 theoreticalKinetic += pseudoParticle[7].GetKineticEnergy() / GeV;
1931 }
1932 }
1933
1934
1935 } else {
1936 theoreticalKinetic -= (currentParticle.GetMass() / GeV + targetParticle.GetMass() / GeV);
1937 for (i = 0; i < vecLen; ++i)
1938 theoreticalKinetic -= vec[i]->GetMass() / GeV;
1939 }
1940 G4double simulatedKinetic = currentParticle.GetKineticEnergy() / GeV + targetParticle.GetKineticEnergy() / GeV;
1941 for (i = 0; i < vecLen; ++i)
1942 simulatedKinetic += vec[i]->GetKineticEnergy() / GeV;
1943
1944
1945
1946
1947
1948 if (simulatedKinetic != 0.0) {
1949 wgt = (theoreticalKinetic) / simulatedKinetic;
1950 currentParticle.SetKineticEnergy(wgt * currentParticle.GetKineticEnergy());
1951 pp = currentParticle.GetTotalMomentum() / MeV;
1952 pp1 = currentParticle.GetMomentum().mag() / MeV;
1953 if (pp1 < 0.001 * MeV) {
1954 rthnve = pi * G4UniformRand();
1955 phinve = twopi * G4UniformRand();
1956 currentParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1957 pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1958 pp * std::cos(rthnve) * MeV);
1959 } else
1960 currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
1961
1962 targetParticle.SetKineticEnergy(wgt * targetParticle.GetKineticEnergy());
1963 pp = targetParticle.GetTotalMomentum() / MeV;
1964 pp1 = targetParticle.GetMomentum().mag() / MeV;
1965 if (pp1 < 0.001 * MeV) {
1966 rthnve = pi * G4UniformRand();
1967 phinve = twopi * G4UniformRand();
1968 targetParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1969 pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1970 pp * std::cos(rthnve) * MeV);
1971 } else
1972 targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
1973
1974 for (i = 0; i < vecLen; ++i) {
1975 vec[i]->SetKineticEnergy(wgt * vec[i]->GetKineticEnergy());
1976 pp = vec[i]->GetTotalMomentum() / MeV;
1977 pp1 = vec[i]->GetMomentum().mag() / MeV;
1978 if (pp1 < 0.001) {
1979 rthnve = pi * G4UniformRand();
1980 phinve = twopi * G4UniformRand();
1981 vec[i]->SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
1982 pp * std::sin(rthnve) * std::sin(phinve) * MeV,
1983 pp * std::cos(rthnve) * MeV);
1984 } else
1985 vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp / pp1));
1986 }
1987 }
1988
1989 Rotate(numberofFinalStateNucleons,
1990 pseudoParticle[4].GetMomentum(),
1991 modifiedOriginal,
1992 originalIncident,
1993 targetNucleus,
1994 currentParticle,
1995 targetParticle,
1996 vec,
1997 vecLen);
1998
1999
2000
2001
2002
2003 if (atomicWeight >= 1.5) {
2004
2005
2006
2007
2008
2009 G4double epnb, edta;
2010 G4int npnb = 0;
2011 G4int ndta = 0;
2012
2013 epnb = targetNucleus.GetPNBlackTrackEnergy();
2014 edta = targetNucleus.GetDTABlackTrackEnergy();
2015 const G4double pnCutOff = 0.001;
2016 const G4double dtaCutOff = 0.001;
2017 const G4double kineticMinimum = 1.e-6;
2018 const G4double kineticFactor = -0.005;
2019
2020 G4double sprob = 0.0;
2021 const G4double ekIncident = originalIncident->GetKineticEnergy() / GeV;
2022 if (ekIncident >= 5.0)
2023 sprob = std::min(1.0, 0.6 * std::log(ekIncident - 4.0));
2024
2025 if (epnb >= pnCutOff) {
2026 npnb = Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * epnb / (epnb + edta));
2027 if (numberofFinalStateNucleons + npnb > atomicWeight)
2028 npnb = G4int(atomicWeight - numberofFinalStateNucleons);
2029 npnb = std::min(npnb, 127 - vecLen);
2030 }
2031 if (edta >= dtaCutOff) {
2032 ndta = Poisson((1.5 + 1.25 * numberofFinalStateNucleons) * edta / (epnb + edta));
2033 ndta = std::min(ndta, 127 - vecLen);
2034 }
2035 G4double spall = numberofFinalStateNucleons;
2036
2037
2038 AddBlackTrackParticles(epnb,
2039 npnb,
2040 edta,
2041 ndta,
2042 sprob,
2043 kineticMinimum,
2044 kineticFactor,
2045 modifiedOriginal,
2046 spall,
2047 targetNucleus,
2048 vec,
2049 vecLen);
2050
2051
2052 }
2053
2054
2055
2056
2057
2058 if ((atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2))
2059 currentParticle.SetTOF(1.0 - 500.0 * std::exp(-ekOriginal / 0.04) * std::log(G4UniformRand()));
2060 else
2061 currentParticle.SetTOF(1.0);
2062
2063
2064 return true;
2065 }
2066
2067 void FullModelReactionDynamics::TwoBody(G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
2068 G4int &vecLen,
2069 G4ReactionProduct &modifiedOriginal,
2070 const G4DynamicParticle * ,
2071 G4ReactionProduct ¤tParticle,
2072 G4ReactionProduct &targetParticle,
2073 const G4Nucleus &targetNucleus,
2074 G4bool & ) {
2075
2076
2077
2078
2079
2080
2081
2082
2083
2084
2085 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
2086 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
2087 G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
2088 G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
2089 G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
2090 G4ParticleDefinition *aKaonZeroS = G4KaonZeroShort::KaonZeroShort();
2091 G4ParticleDefinition *aKaonZeroL = G4KaonZeroLong::KaonZeroLong();
2092
2093
2094 static const G4double expxu = 82.;
2095 static const G4double expxl = -expxu;
2096
2097 const G4double ekOriginal = modifiedOriginal.GetKineticEnergy() / GeV;
2098 const G4double pOriginal = modifiedOriginal.GetMomentum().mag() / GeV;
2099 G4double currentMass = currentParticle.GetMass() / GeV;
2100 G4double targetMass = targetParticle.GetMass() / GeV;
2101 const G4double atomicWeight = targetNucleus.GetN_asInt();
2102
2103 G4double etCurrent = currentParticle.GetTotalEnergy() / GeV;
2104 G4double pCurrent = currentParticle.GetTotalMomentum() / GeV;
2105
2106 G4double cmEnergy =
2107 std::sqrt(currentMass * currentMass + targetMass * targetMass + 2.0 * targetMass * etCurrent);
2108
2109
2110
2111
2112
2113
2114
2115 if ((pCurrent < 0.1) || (cmEnergy < 0.01))
2116 {
2117 targetParticle.SetMass(0.0);
2118 } else {
2119
2120
2121
2122
2123
2124
2125
2126
2127
2128
2129
2130
2131
2132
2133
2134
2135
2136
2137
2138
2139
2140
2141
2142
2143
2144
2145
2146
2147
2148 G4double pf = cmEnergy * cmEnergy + targetMass * targetMass - currentMass * currentMass;
2149 pf = pf * pf - 4 * cmEnergy * cmEnergy * targetMass * targetMass;
2150
2151
2152 if (pf <= 0.)
2153 {
2154 for (G4int i = 0; i < vecLen; i++)
2155 delete vec[i];
2156 vecLen = 0;
2157 throw G4HadronicException(__FILE__, __LINE__, "FullModelReactionDynamics::TwoBody: pf is too small ");
2158 }
2159
2160 pf = std::sqrt(pf) / (2.0 * cmEnergy);
2161
2162
2163
2164 G4ReactionProduct pseudoParticle[3];
2165 pseudoParticle[0].SetMass(currentMass * GeV);
2166 pseudoParticle[0].SetTotalEnergy(etCurrent * GeV);
2167 pseudoParticle[0].SetMomentum(0.0, 0.0, pCurrent * GeV);
2168
2169 pseudoParticle[1].SetMomentum(0.0, 0.0, 0.0);
2170 pseudoParticle[1].SetMass(targetMass * GeV);
2171 pseudoParticle[1].SetKineticEnergy(0.0);
2172
2173
2174
2175 pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
2176 pseudoParticle[0].Lorentz(pseudoParticle[0], pseudoParticle[2]);
2177 pseudoParticle[1].Lorentz(pseudoParticle[1], pseudoParticle[2]);
2178
2179
2180
2181 currentParticle.SetTotalEnergy(std::sqrt(pf * pf + currentMass * currentMass) * GeV);
2182 targetParticle.SetTotalEnergy(std::sqrt(pf * pf + targetMass * targetMass) * GeV);
2183
2184
2185
2186 const G4double cb = 0.01;
2187 const G4double b1 = 4.225;
2188 const G4double b2 = 1.795;
2189
2190
2191
2192 G4double b = std::max(cb, b1 + b2 * std::log(pOriginal));
2193
2194 G4double btrang = b * 4.0 * pf * pseudoParticle[0].GetMomentum().mag() / GeV;
2195
2196 G4double exindt = -1.0;
2197 exindt += std::exp(std::max(-btrang, expxl));
2198
2199
2200
2201 G4double ctet = 1.0 + 2 * std::log(1.0 + G4UniformRand() * exindt) / btrang;
2202 if (std::fabs(ctet) > 1.0)
2203 ctet > 0.0 ? ctet = 1.0 : ctet = -1.0;
2204 G4double stet = std::sqrt((1.0 - ctet) * (1.0 + ctet));
2205 G4double phi = twopi * G4UniformRand();
2206
2207
2208
2209 if (targetParticle.GetDefinition() == aKaonMinus || targetParticle.GetDefinition() == aKaonZeroL ||
2210 targetParticle.GetDefinition() == aKaonZeroS || targetParticle.GetDefinition() == aKaonPlus ||
2211 targetParticle.GetDefinition() == aPiMinus || targetParticle.GetDefinition() == aPiZero ||
2212 targetParticle.GetDefinition() == aPiPlus) {
2213 currentParticle.SetMomentum(-pf * stet * std::sin(phi) * GeV, -pf * stet * std::cos(phi) * GeV, -pf * ctet * GeV);
2214 } else {
2215 currentParticle.SetMomentum(pf * stet * std::sin(phi) * GeV, pf * stet * std::cos(phi) * GeV, pf * ctet * GeV);
2216 }
2217 targetParticle.SetMomentum(currentParticle.GetMomentum() * (-1.0));
2218
2219
2220
2221 currentParticle.Lorentz(currentParticle, pseudoParticle[1]);
2222 targetParticle.Lorentz(targetParticle, pseudoParticle[1]);
2223
2224
2225
2226
2227
2228
2229
2230
2231
2232 Defs1(modifiedOriginal, currentParticle, targetParticle, vec, vecLen);
2233
2234 G4double pp, pp1, ekin;
2235 if (atomicWeight >= 1.5) {
2236 const G4double cfa = 0.025 * ((atomicWeight - 1.) / 120.) * std::exp(-(atomicWeight - 1.) / 120.);
2237 pp1 = currentParticle.GetMomentum().mag() / MeV;
2238 if (pp1 >= 1.0) {
2239 ekin = currentParticle.GetKineticEnergy() / MeV - cfa * (1.0 + 0.5 * normal()) * GeV;
2240 ekin = std::max(0.0001 * GeV, ekin);
2241 currentParticle.SetKineticEnergy(ekin * MeV);
2242 pp = currentParticle.GetTotalMomentum() / MeV;
2243 currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
2244 }
2245 pp1 = targetParticle.GetMomentum().mag() / MeV;
2246 if (pp1 >= 1.0) {
2247 ekin = targetParticle.GetKineticEnergy() / MeV - cfa * (1.0 + normal() / 2.) * GeV;
2248 ekin = std::max(0.0001 * GeV, ekin);
2249 targetParticle.SetKineticEnergy(ekin * MeV);
2250 pp = targetParticle.GetTotalMomentum() / MeV;
2251 targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
2252 }
2253 }
2254 }
2255
2256 if (atomicWeight >= 1.5) {
2257
2258
2259
2260
2261
2262
2263
2264
2265
2266
2267 G4double epnb, edta;
2268 G4int npnb = 0, ndta = 0;
2269
2270 epnb = targetNucleus.GetPNBlackTrackEnergy();
2271 edta = targetNucleus.GetDTABlackTrackEnergy();
2272 const G4double pnCutOff = 0.0001;
2273 const G4double dtaCutOff = 0.0001;
2274 const G4double kineticMinimum = 0.0001;
2275 const G4double kineticFactor = -0.010;
2276 G4double sprob = 0.0;
2277 if (epnb >= pnCutOff) {
2278 npnb = Poisson(epnb / 0.02);
2279
2280
2281
2282
2283
2284 if (npnb > atomicWeight)
2285 npnb = G4int(atomicWeight);
2286 if ((epnb > pnCutOff) && (npnb <= 0))
2287 npnb = 1;
2288 npnb = std::min(npnb, 127 - vecLen);
2289 }
2290 if (edta >= dtaCutOff) {
2291 ndta = G4int(2.0 * std::log(atomicWeight));
2292 ndta = std::min(ndta, 127 - vecLen);
2293 }
2294 G4double spall = 0.0;
2295
2296
2297
2298
2299
2300
2301
2302
2303
2304
2305
2306
2307
2308
2309
2310
2311
2312
2313 AddBlackTrackParticles(epnb,
2314 npnb,
2315 edta,
2316 ndta,
2317 sprob,
2318 kineticMinimum,
2319 kineticFactor,
2320 modifiedOriginal,
2321 spall,
2322 targetNucleus,
2323 vec,
2324 vecLen);
2325
2326
2327 }
2328
2329
2330
2331 if ((atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2))
2332 currentParticle.SetTOF(1.0 - 500.0 * std::exp(-ekOriginal / 0.04) * std::log(G4UniformRand()));
2333 else
2334 currentParticle.SetTOF(1.0);
2335 return;
2336 }
2337
2338 G4double FullModelReactionDynamics::GenerateNBodyEvent(const G4double totalEnergy,
2339 const G4bool constantCrossSection,
2340 G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
2341 G4int &vecLen) {
2342
2343
2344
2345
2346 G4int i;
2347 const G4double expxu = 82.;
2348 const G4double expxl = -expxu;
2349 if (vecLen < 2) {
2350 G4cerr << "*** Error in FullModelReactionDynamics::GenerateNBodyEvent" << G4endl;
2351 G4cerr << " number of particles < 2" << G4endl;
2352 G4cerr << "totalEnergy = " << totalEnergy << "MeV, vecLen = " << vecLen << G4endl;
2353 return -1.0;
2354 }
2355 G4double mass[18];
2356 G4double energy[18];
2357 G4double pcm[3][18];
2358 G4double totalMass = 0.0;
2359
2360 G4double sm[18];
2361
2362 for (i = 0; i < vecLen; ++i) {
2363 mass[i] = vec[i]->GetMass() / GeV;
2364
2365
2366 vec[i]->SetMomentum(0.0, 0.0, 0.0);
2367 pcm[0][i] = 0.0;
2368 pcm[1][i] = 0.0;
2369 pcm[2][i] = 0.0;
2370 energy[i] = mass[i];
2371 totalMass += mass[i];
2372 sm[i] = totalMass;
2373 }
2374 G4double totalE = totalEnergy / GeV;
2375 if (totalMass > totalE) {
2376 return -1.0;
2377 }
2378 G4double kineticEnergy = totalE - totalMass;
2379 G4double emm[18];
2380
2381 emm[0] = mass[0];
2382 emm[vecLen - 1] = totalE;
2383 if (vecLen > 2)
2384 {
2385 G4double ran[18];
2386 for (i = 0; i < vecLen; ++i)
2387 ran[i] = G4UniformRand();
2388 for (i = 0; i < vecLen - 2; ++i) {
2389 for (G4int j = vecLen - 2; j > i; --j) {
2390 if (ran[i] > ran[j]) {
2391 G4double temp = ran[i];
2392 ran[i] = ran[j];
2393 ran[j] = temp;
2394 }
2395 }
2396 }
2397 for (i = 1; i < vecLen - 1; ++i)
2398 emm[i] = ran[i - 1] * kineticEnergy + sm[i];
2399 }
2400
2401 G4bool lzero = true;
2402 G4double wtmax = 0.0;
2403 if (constantCrossSection)
2404 {
2405 G4double emmax = kineticEnergy + mass[0];
2406 G4double emmin = 0.0;
2407 for (i = 1; i < vecLen; ++i) {
2408 emmin += mass[i - 1];
2409 emmax += mass[i];
2410 G4double wtfc = 0.0;
2411 if (emmax * emmax > 0.0) {
2412 G4double arg = emmax * emmax +
2413 (emmin * emmin - mass[i] * mass[i]) * (emmin * emmin - mass[i] * mass[i]) / (emmax * emmax) -
2414 2.0 * (emmin * emmin + mass[i] * mass[i]);
2415 if (arg > 0.0)
2416 wtfc = 0.5 * std::sqrt(arg);
2417 }
2418 if (wtfc == 0.0) {
2419 lzero = false;
2420 break;
2421 }
2422 wtmax += std::log(wtfc);
2423 }
2424 if (lzero)
2425 wtmax = -wtmax;
2426 else
2427 wtmax = expxu;
2428 } else {
2429
2430 const G4double ffq[18] = {0.,
2431 3.141592,
2432 19.73921,
2433 62.01255,
2434 129.8788,
2435 204.0131,
2436 256.3704,
2437 268.4705,
2438 240.9780,
2439 189.2637,
2440 132.1308,
2441 83.0202,
2442 47.4210,
2443 24.8295,
2444 12.0006,
2445 5.3858,
2446 2.2560,
2447 0.8859};
2448 wtmax = std::log(std::pow(kineticEnergy, vecLen - 2) * ffq[vecLen - 1] / totalE);
2449 }
2450 lzero = true;
2451 G4double pd[50];
2452
2453 for (i = 0; i < vecLen - 1; ++i) {
2454 pd[i] = 0.0;
2455 if (emm[i + 1] * emm[i + 1] > 0.0) {
2456 G4double arg = emm[i + 1] * emm[i + 1] +
2457 (emm[i] * emm[i] - mass[i + 1] * mass[i + 1]) * (emm[i] * emm[i] - mass[i + 1] * mass[i + 1]) /
2458 (emm[i + 1] * emm[i + 1]) -
2459 2.0 * (emm[i] * emm[i] + mass[i + 1] * mass[i + 1]);
2460 if (arg > 0.0)
2461 pd[i] = 0.5 * std::sqrt(arg);
2462 }
2463 if (pd[i] <= 0.0)
2464 lzero = false;
2465 else
2466 wtmax += std::log(pd[i]);
2467 }
2468 G4double weight = 0.0;
2469 if (lzero)
2470 weight = std::exp(std::max(std::min(wtmax, expxu), expxl));
2471
2472 G4double bang, cb, sb, s0, s1, s2, c, s, esys, a, b, gama, beta;
2473 pcm[0][0] = 0.0;
2474 pcm[1][0] = pd[0];
2475 pcm[2][0] = 0.0;
2476 for (i = 1; i < vecLen; ++i) {
2477 pcm[0][i] = 0.0;
2478 pcm[1][i] = -pd[i - 1];
2479 pcm[2][i] = 0.0;
2480 bang = twopi * G4UniformRand();
2481 cb = std::cos(bang);
2482 sb = std::sin(bang);
2483 c = 2.0 * G4UniformRand() - 1.0;
2484 s = std::sqrt(std::fabs(1.0 - c * c));
2485 if (i < vecLen - 1) {
2486 esys = std::sqrt(pd[i] * pd[i] + emm[i] * emm[i]);
2487 beta = pd[i] / esys;
2488 gama = esys / emm[i];
2489 for (G4int j = 0; j <= i; ++j) {
2490 s0 = pcm[0][j];
2491 s1 = pcm[1][j];
2492 s2 = pcm[2][j];
2493 energy[j] = std::sqrt(s0 * s0 + s1 * s1 + s2 * s2 + mass[j] * mass[j]);
2494 a = s0 * c - s1 * s;
2495 pcm[1][j] = s0 * s + s1 * c;
2496 b = pcm[2][j];
2497 pcm[0][j] = a * cb - b * sb;
2498 pcm[2][j] = a * sb + b * cb;
2499 pcm[1][j] = gama * (pcm[1][j] + beta * energy[j]);
2500 }
2501 } else {
2502 for (G4int j = 0; j <= i; ++j) {
2503 s0 = pcm[0][j];
2504 s1 = pcm[1][j];
2505 s2 = pcm[2][j];
2506 energy[j] = std::sqrt(s0 * s0 + s1 * s1 + s2 * s2 + mass[j] * mass[j]);
2507 a = s0 * c - s1 * s;
2508 pcm[1][j] = s0 * s + s1 * c;
2509 b = pcm[2][j];
2510 pcm[0][j] = a * cb - b * sb;
2511 pcm[2][j] = a * sb + b * cb;
2512 }
2513 }
2514 }
2515 for (i = 0; i < vecLen; ++i) {
2516 vec[i]->SetMomentum(pcm[0][i] * GeV, pcm[1][i] * GeV, pcm[2][i] * GeV);
2517 vec[i]->SetTotalEnergy(energy[i] * GeV);
2518 }
2519
2520 return weight;
2521 }
2522
2523 G4double FullModelReactionDynamics::normal() {
2524 G4double ran = -6.0;
2525 for (G4int i = 0; i < 12; ++i)
2526 ran += G4UniformRand();
2527 return ran;
2528 }
2529
2530 G4int FullModelReactionDynamics::Poisson(G4double x)
2531 {
2532 G4int iran;
2533 G4double ran;
2534
2535 if (x > 9.9)
2536 iran = static_cast<G4int>(std::max(0.0, x + normal() * std::sqrt(x)));
2537 else {
2538 G4int mm = G4int(5.0 * x);
2539 if (mm <= 0)
2540 {
2541 G4double p1 = x * std::exp(-x);
2542 G4double p2 = x * p1 / 2.0;
2543 G4double p3 = x * p2 / 3.0;
2544 ran = G4UniformRand();
2545 if (ran < p3)
2546 iran = 3;
2547 else if (ran < p2)
2548 iran = 2;
2549 else if (ran < p1)
2550 iran = 1;
2551 else
2552 iran = 0;
2553 } else {
2554 iran = 0;
2555 G4double r = std::exp(-x);
2556 ran = G4UniformRand();
2557 if (ran > r) {
2558 G4double rrr;
2559 G4double rr = r;
2560 for (G4int i = 1; i <= mm; ++i) {
2561 iran++;
2562 if (i > 5)
2563 rrr = std::exp(i * std::log(x) - (i + 0.5) * std::log((G4double)i) + i - 0.9189385);
2564 else
2565 rrr = std::pow(x, i) / Factorial(i);
2566 rr += r * rrr;
2567 if (ran <= rr)
2568 break;
2569 }
2570 }
2571 }
2572 }
2573 return iran;
2574 }
2575
2576 G4int FullModelReactionDynamics::Factorial(G4int n) {
2577 G4int m = std::min(n, 10);
2578 G4int result = 1;
2579 if (m <= 1)
2580 return result;
2581 for (G4int i = 2; i <= m; ++i)
2582 result *= i;
2583 return result;
2584 }
2585
2586 void FullModelReactionDynamics::Defs1(const G4ReactionProduct &modifiedOriginal,
2587 G4ReactionProduct ¤tParticle,
2588 G4ReactionProduct &targetParticle,
2589 G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
2590 G4int &vecLen) {
2591 const G4double pjx = modifiedOriginal.GetMomentum().x() / MeV;
2592 const G4double pjy = modifiedOriginal.GetMomentum().y() / MeV;
2593 const G4double pjz = modifiedOriginal.GetMomentum().z() / MeV;
2594 const G4double p = modifiedOriginal.GetMomentum().mag() / MeV;
2595 if (pjx * pjx + pjy * pjy > 0.0) {
2596 G4double cost, sint, ph, cosp, sinp, pix, piy, piz;
2597 cost = pjz / p;
2598 sint = 0.5 * (std::sqrt(std::abs((1.0 - cost) * (1.0 + cost))) + std::sqrt(pjx * pjx + pjy * pjy) / p);
2599 if (pjy < 0.0)
2600 ph = 3 * halfpi;
2601 else
2602 ph = halfpi;
2603 if (std::abs(pjx) > 0.001 * MeV)
2604 ph = std::atan2(pjy, pjx);
2605 cosp = std::cos(ph);
2606 sinp = std::sin(ph);
2607 pix = currentParticle.GetMomentum().x() / MeV;
2608 piy = currentParticle.GetMomentum().y() / MeV;
2609 piz = currentParticle.GetMomentum().z() / MeV;
2610 currentParticle.SetMomentum(cost * cosp * pix * MeV - sinp * piy + sint * cosp * piz * MeV,
2611 cost * sinp * pix * MeV + cosp * piy + sint * sinp * piz * MeV,
2612 -sint * pix * MeV + cost * piz * MeV);
2613 pix = targetParticle.GetMomentum().x() / MeV;
2614 piy = targetParticle.GetMomentum().y() / MeV;
2615 piz = targetParticle.GetMomentum().z() / MeV;
2616 targetParticle.SetMomentum(cost * cosp * pix * MeV - sinp * piy + sint * cosp * piz * MeV,
2617 cost * sinp * pix * MeV + cosp * piy + sint * sinp * piz * MeV,
2618 -sint * pix * MeV + cost * piz * MeV);
2619 for (G4int i = 0; i < vecLen; ++i) {
2620 pix = vec[i]->GetMomentum().x() / MeV;
2621 piy = vec[i]->GetMomentum().y() / MeV;
2622 piz = vec[i]->GetMomentum().z() / MeV;
2623 vec[i]->SetMomentum(cost * cosp * pix * MeV - sinp * piy + sint * cosp * piz * MeV,
2624 cost * sinp * pix * MeV + cosp * piy + sint * sinp * piz * MeV,
2625 -sint * pix * MeV + cost * piz * MeV);
2626 }
2627 } else {
2628 if (pjz < 0.0) {
2629 currentParticle.SetMomentum(-currentParticle.GetMomentum().z());
2630 targetParticle.SetMomentum(-targetParticle.GetMomentum().z());
2631 for (G4int i = 0; i < vecLen; ++i)
2632 vec[i]->SetMomentum(-vec[i]->GetMomentum().z());
2633 }
2634 }
2635 }
2636
2637 void FullModelReactionDynamics::Rotate(
2638 const G4double numberofFinalStateNucleons,
2639 const G4ThreeVector &temp,
2640 const G4ReactionProduct &modifiedOriginal,
2641 const G4HadProjectile *originalIncident,
2642 const G4Nucleus &targetNucleus,
2643 G4ReactionProduct ¤tParticle,
2644 G4ReactionProduct &targetParticle,
2645 G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
2646 G4int &vecLen) {
2647
2648
2649
2650
2651
2652 const G4double atomicWeight = targetNucleus.GetN_asInt();
2653 const G4double logWeight = std::log(atomicWeight);
2654
2655 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
2656 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
2657 G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
2658
2659 G4int i;
2660 G4ThreeVector pseudoParticle[4];
2661 for (i = 0; i < 4; ++i)
2662 pseudoParticle[i].set(0, 0, 0);
2663 pseudoParticle[0] = currentParticle.GetMomentum() + targetParticle.GetMomentum();
2664 for (i = 0; i < vecLen; ++i)
2665 pseudoParticle[0] = pseudoParticle[0] + (vec[i]->GetMomentum());
2666
2667
2668
2669 G4float pp, pp1;
2670 G4double alekw, p, rthnve, phinve;
2671 G4double r1, r2, a1, ran1, ran2, xxh, exh, pxTemp, pyTemp, pzTemp;
2672
2673 r1 = twopi * G4UniformRand();
2674 r2 = G4UniformRand();
2675 a1 = std::sqrt(-2.0 * std::log(r2));
2676 ran1 = a1 * std::sin(r1) * 0.020 * numberofFinalStateNucleons * GeV;
2677 ran2 = a1 * std::cos(r1) * 0.020 * numberofFinalStateNucleons * GeV;
2678 G4ThreeVector fermi(ran1, ran2, 0);
2679
2680 pseudoParticle[0] = pseudoParticle[0] + fermi;
2681 pseudoParticle[2] = temp;
2682 pseudoParticle[3] = pseudoParticle[0];
2683
2684 pseudoParticle[1] = pseudoParticle[2].cross(pseudoParticle[3]);
2685 G4double rotation = 2. * pi * G4UniformRand();
2686 pseudoParticle[1] = pseudoParticle[1].rotate(rotation, pseudoParticle[3]);
2687 pseudoParticle[2] = pseudoParticle[3].cross(pseudoParticle[1]);
2688 for (G4int ii = 1; ii <= 3; ii++) {
2689 p = pseudoParticle[ii].mag();
2690 if (p == 0.0)
2691 pseudoParticle[ii] = G4ThreeVector(0.0, 0.0, 0.0);
2692 else
2693 pseudoParticle[ii] = pseudoParticle[ii] * (1. / p);
2694 }
2695
2696 pxTemp = pseudoParticle[1].dot(currentParticle.GetMomentum());
2697 pyTemp = pseudoParticle[2].dot(currentParticle.GetMomentum());
2698 pzTemp = pseudoParticle[3].dot(currentParticle.GetMomentum());
2699 currentParticle.SetMomentum(pxTemp, pyTemp, pzTemp);
2700
2701 pxTemp = pseudoParticle[1].dot(targetParticle.GetMomentum());
2702 pyTemp = pseudoParticle[2].dot(targetParticle.GetMomentum());
2703 pzTemp = pseudoParticle[3].dot(targetParticle.GetMomentum());
2704 targetParticle.SetMomentum(pxTemp, pyTemp, pzTemp);
2705
2706 for (i = 0; i < vecLen; ++i) {
2707 pxTemp = pseudoParticle[1].dot(vec[i]->GetMomentum());
2708 pyTemp = pseudoParticle[2].dot(vec[i]->GetMomentum());
2709 pzTemp = pseudoParticle[3].dot(vec[i]->GetMomentum());
2710 vec[i]->SetMomentum(pxTemp, pyTemp, pzTemp);
2711 }
2712
2713
2714
2715
2716 Defs1(modifiedOriginal, currentParticle, targetParticle, vec, vecLen);
2717 G4double ekin;
2718 G4double dekin = 0.0;
2719 G4double ek1 = 0.0;
2720 G4int npions = 0;
2721 if (atomicWeight >= 1.5)
2722 {
2723
2724
2725 const G4double alem[] = {1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00};
2726 const G4double val0[] = {0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65};
2727 alekw = std::log(originalIncident->GetKineticEnergy() / GeV);
2728 exh = 1.0;
2729 if (alekw > alem[0])
2730 {
2731 exh = val0[6];
2732 for (G4int j = 1; j < 7; ++j) {
2733 if (alekw < alem[j])
2734 {
2735 G4double rcnve = (val0[j] - val0[j - 1]) / (alem[j] - alem[j - 1]);
2736 exh = rcnve * alekw + val0[j - 1] - rcnve * alem[j - 1];
2737 break;
2738 }
2739 }
2740 exh = 1.0 - exh;
2741 }
2742 const G4double cfa = 0.025 * ((atomicWeight - 1.) / 120.) * std::exp(-(atomicWeight - 1.) / 120.);
2743 ekin = currentParticle.GetKineticEnergy() / GeV - cfa * (1 + normal() / 2.0);
2744 ekin = std::max(1.0e-6, ekin);
2745 xxh = 1.0;
2746 dekin += ekin * (1.0 - xxh);
2747 ekin *= xxh;
2748 currentParticle.SetKineticEnergy(ekin * GeV);
2749 pp = currentParticle.GetTotalMomentum() / MeV;
2750 pp1 = currentParticle.GetMomentum().mag() / MeV;
2751 if (pp1 < 0.001 * MeV) {
2752 rthnve = pi * G4UniformRand();
2753 phinve = twopi * G4UniformRand();
2754 currentParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2755 pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2756 pp * std::cos(rthnve) * MeV);
2757 } else
2758 currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
2759 ekin = targetParticle.GetKineticEnergy() / GeV - cfa * (1 + normal() / 2.0);
2760 ekin = std::max(1.0e-6, ekin);
2761 xxh = 1.0;
2762 if (((modifiedOriginal.GetDefinition() == aPiPlus) || (modifiedOriginal.GetDefinition() == aPiMinus)) &&
2763 (targetParticle.GetDefinition() == aPiZero) && (G4UniformRand() < logWeight))
2764 xxh = exh;
2765 dekin += ekin * (1.0 - xxh);
2766 ekin *= xxh;
2767 if ((targetParticle.GetDefinition() == aPiPlus) || (targetParticle.GetDefinition() == aPiZero) ||
2768 (targetParticle.GetDefinition() == aPiMinus)) {
2769 ++npions;
2770 ek1 += ekin;
2771 }
2772 targetParticle.SetKineticEnergy(ekin * GeV);
2773 pp = targetParticle.GetTotalMomentum() / MeV;
2774 pp1 = targetParticle.GetMomentum().mag() / MeV;
2775 if (pp1 < 0.001 * MeV) {
2776 rthnve = pi * G4UniformRand();
2777 phinve = twopi * G4UniformRand();
2778 targetParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2779 pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2780 pp * std::cos(rthnve) * MeV);
2781 } else
2782 targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp / pp1));
2783 for (i = 0; i < vecLen; ++i) {
2784 ekin = vec[i]->GetKineticEnergy() / GeV - cfa * (1 + normal() / 2.0);
2785 ekin = std::max(1.0e-6, ekin);
2786 xxh = 1.0;
2787 if (((modifiedOriginal.GetDefinition() == aPiPlus) || (modifiedOriginal.GetDefinition() == aPiMinus)) &&
2788 (vec[i]->GetDefinition() == aPiZero) && (G4UniformRand() < logWeight))
2789 xxh = exh;
2790 dekin += ekin * (1.0 - xxh);
2791 ekin *= xxh;
2792 if ((vec[i]->GetDefinition() == aPiPlus) || (vec[i]->GetDefinition() == aPiZero) ||
2793 (vec[i]->GetDefinition() == aPiMinus)) {
2794 ++npions;
2795 ek1 += ekin;
2796 }
2797 vec[i]->SetKineticEnergy(ekin * GeV);
2798 pp = vec[i]->GetTotalMomentum() / MeV;
2799 pp1 = vec[i]->GetMomentum().mag() / MeV;
2800 if (pp1 < 0.001 * MeV) {
2801 rthnve = pi * G4UniformRand();
2802 phinve = twopi * G4UniformRand();
2803 vec[i]->SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2804 pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2805 pp * std::cos(rthnve) * MeV);
2806 } else
2807 vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp / pp1));
2808 }
2809 }
2810 if ((ek1 != 0.0) && (npions > 0)) {
2811 dekin = 1.0 + dekin / ek1;
2812
2813
2814
2815 if ((currentParticle.GetDefinition() == aPiPlus) || (currentParticle.GetDefinition() == aPiZero) ||
2816 (currentParticle.GetDefinition() == aPiMinus)) {
2817 currentParticle.SetKineticEnergy(std::max(0.001 * MeV, dekin * currentParticle.GetKineticEnergy()));
2818 pp = currentParticle.GetTotalMomentum() / MeV;
2819 pp1 = currentParticle.GetMomentum().mag() / MeV;
2820 if (pp1 < 0.001) {
2821 rthnve = pi * G4UniformRand();
2822 phinve = twopi * G4UniformRand();
2823 currentParticle.SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2824 pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2825 pp * std::cos(rthnve) * MeV);
2826 } else
2827 currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp / pp1));
2828 }
2829 for (i = 0; i < vecLen; ++i) {
2830 if ((vec[i]->GetDefinition() == aPiPlus) || (vec[i]->GetDefinition() == aPiZero) ||
2831 (vec[i]->GetDefinition() == aPiMinus)) {
2832 vec[i]->SetKineticEnergy(std::max(0.001 * MeV, dekin * vec[i]->GetKineticEnergy()));
2833 pp = vec[i]->GetTotalMomentum() / MeV;
2834 pp1 = vec[i]->GetMomentum().mag() / MeV;
2835 if (pp1 < 0.001) {
2836 rthnve = pi * G4UniformRand();
2837 phinve = twopi * G4UniformRand();
2838 vec[i]->SetMomentum(pp * std::sin(rthnve) * std::cos(phinve) * MeV,
2839 pp * std::sin(rthnve) * std::sin(phinve) * MeV,
2840 pp * std::cos(rthnve) * MeV);
2841 } else
2842 vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp / pp1));
2843 }
2844 }
2845 }
2846 }
2847
2848 void FullModelReactionDynamics::AddBlackTrackParticles(const G4double epnb,
2849 const G4int npnb,
2850 const G4double edta,
2851 const G4int ndta,
2852 const G4double sprob,
2853 const G4double kineticMinimum,
2854 const G4double kineticFactor,
2855 const G4ReactionProduct &modifiedOriginal,
2856 G4double spall,
2857 const G4Nucleus &targetNucleus,
2858 G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
2859 G4int &vecLen) {
2860
2861
2862
2863
2864
2865
2866
2867
2868 G4ParticleDefinition *aProton = G4Proton::Proton();
2869 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
2870 G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
2871 G4ParticleDefinition *aTriton = G4Triton::Triton();
2872 G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
2873
2874 const G4double ekOriginal = modifiedOriginal.GetKineticEnergy() / MeV;
2875 G4double atomicWeight = targetNucleus.GetN_asInt();
2876 G4double atomicNumber = targetNucleus.GetZ_asInt();
2877
2878 const G4double ika1 = 3.6;
2879 const G4double ika2 = 35.56;
2880 const G4double ika3 = 6.45;
2881 const G4double sp1 = 1.066;
2882
2883 G4int i;
2884 G4double pp;
2885
2886
2887 G4double cfa = 0.025 * ((atomicWeight - 1.0) / 120.0) * std::exp(-(atomicWeight - 1.0) / 120.0);
2888 if (npnb > 0)
2889 {
2890 G4double backwardKinetic = 0.0;
2891 G4int local_npnb = npnb;
2892 for (i = 0; i < npnb; ++i)
2893 if (G4UniformRand() < sprob)
2894 local_npnb--;
2895 G4double ekin = epnb / local_npnb;
2896
2897 for (i = 0; i < local_npnb; ++i) {
2898 G4ReactionProduct *p1 = new G4ReactionProduct();
2899 if (backwardKinetic > epnb) {
2900 delete p1;
2901 break;
2902 }
2903 G4double ran = G4UniformRand();
2904 G4double kinetic = -ekin * std::log(ran) - cfa * (1.0 + 0.5 * normal());
2905 if (kinetic < 0.0)
2906 kinetic = -0.010 * std::log(ran);
2907 backwardKinetic += kinetic;
2908 if (backwardKinetic > epnb)
2909 kinetic = std::max(kineticMinimum, epnb - (backwardKinetic - kinetic));
2910 if (G4UniformRand() > (1.0 - atomicNumber / atomicWeight))
2911 p1->SetDefinition(aProton);
2912 else
2913 p1->SetDefinition(aNeutron);
2914 vec.SetElement(vecLen, p1);
2915 ++spall;
2916 G4double cost = G4UniformRand() * 2.0 - 1.0;
2917 G4double sint = std::sqrt(std::fabs(1.0 - cost * cost));
2918 G4double phi = twopi * G4UniformRand();
2919 vec[vecLen]->SetNewlyAdded(true);
2920 vec[vecLen]->SetKineticEnergy(kinetic * GeV);
2921
2922 pp = vec[vecLen]->GetTotalMomentum() / MeV;
2923 vec[vecLen]->SetMomentum(pp * sint * std::sin(phi) * MeV, pp * sint * std::cos(phi) * MeV, pp * cost * MeV);
2924 vecLen++;
2925
2926 }
2927 if ((atomicWeight >= 10.0) && (ekOriginal <= 2.0 * GeV)) {
2928 G4double ekw = ekOriginal / GeV;
2929 G4int ika, kk = 0;
2930 if (ekw > 1.0)
2931 ekw *= ekw;
2932 ekw = std::max(0.1, ekw);
2933 ika = G4int(ika1 * std::exp((atomicNumber * atomicNumber / atomicWeight - ika2) / ika3) / ekw);
2934 if (ika > 0) {
2935 for (i = (vecLen - 1); i >= 0; --i) {
2936 if ((vec[i]->GetDefinition() == aProton) && vec[i]->GetNewlyAdded()) {
2937 vec[i]->SetDefinitionAndUpdateE(aNeutron);
2938 if (++kk > ika)
2939 break;
2940 }
2941 }
2942 }
2943 }
2944 }
2945 if (ndta > 0)
2946 {
2947 G4double backwardKinetic = 0.0;
2948 G4int local_ndta = ndta;
2949 for (i = 0; i < ndta; ++i)
2950 if (G4UniformRand() < sprob)
2951 local_ndta--;
2952 G4double ekin = edta / local_ndta;
2953
2954 for (i = 0; i < local_ndta; ++i) {
2955 G4ReactionProduct *p2 = new G4ReactionProduct();
2956 if (backwardKinetic > edta) {
2957 delete p2;
2958 break;
2959 }
2960 G4double ran = G4UniformRand();
2961 G4double kinetic = -ekin * std::log(ran) - cfa * (1. + 0.5 * normal());
2962 if (kinetic < 0.0)
2963 kinetic = kineticFactor * std::log(ran);
2964 backwardKinetic += kinetic;
2965 if (backwardKinetic > edta)
2966 kinetic = edta - (backwardKinetic - kinetic);
2967 if (kinetic < 0.0)
2968 kinetic = kineticMinimum;
2969 G4double cost = 2.0 * G4UniformRand() - 1.0;
2970 G4double sint = std::sqrt(std::max(0.0, (1.0 - cost * cost)));
2971 G4double phi = twopi * G4UniformRand();
2972 ran = G4UniformRand();
2973 if (ran <= 0.60)
2974 p2->SetDefinition(aDeuteron);
2975 else if (ran <= 0.90)
2976 p2->SetDefinition(aTriton);
2977 else
2978 p2->SetDefinition(anAlpha);
2979 spall += p2->GetMass() / GeV * sp1;
2980 if (spall > atomicWeight) {
2981 delete p2;
2982 break;
2983 }
2984 vec.SetElement(vecLen, p2);
2985 vec[vecLen]->SetNewlyAdded(true);
2986 vec[vecLen]->SetKineticEnergy(kinetic * GeV);
2987
2988 pp = vec[vecLen]->GetTotalMomentum() / MeV;
2989 vec[vecLen++]->SetMomentum(pp * sint * std::sin(phi) * MeV, pp * sint * std::cos(phi) * MeV, pp * cost * MeV);
2990
2991 }
2992 }
2993
2994 }
2995
2996 void FullModelReactionDynamics::MomentumCheck(const G4ReactionProduct &modifiedOriginal,
2997 G4ReactionProduct ¤tParticle,
2998 G4ReactionProduct &targetParticle,
2999 G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
3000 G4int &vecLen) {
3001 const G4double pOriginal = modifiedOriginal.GetTotalMomentum() / MeV;
3002 G4double testMomentum = currentParticle.GetMomentum().mag() / MeV;
3003 G4double pMass;
3004 if (testMomentum >= pOriginal) {
3005 pMass = currentParticle.GetMass() / MeV;
3006 currentParticle.SetTotalEnergy(std::sqrt(pMass * pMass + pOriginal * pOriginal) * MeV);
3007 currentParticle.SetMomentum(currentParticle.GetMomentum() * (pOriginal / testMomentum));
3008 }
3009 testMomentum = targetParticle.GetMomentum().mag() / MeV;
3010 if (testMomentum >= pOriginal) {
3011 pMass = targetParticle.GetMass() / MeV;
3012 targetParticle.SetTotalEnergy(std::sqrt(pMass * pMass + pOriginal * pOriginal) * MeV);
3013 targetParticle.SetMomentum(targetParticle.GetMomentum() * (pOriginal / testMomentum));
3014 }
3015 for (G4int i = 0; i < vecLen; ++i) {
3016 testMomentum = vec[i]->GetMomentum().mag() / MeV;
3017 if (testMomentum >= pOriginal) {
3018 pMass = vec[i]->GetMass() / MeV;
3019 vec[i]->SetTotalEnergy(std::sqrt(pMass * pMass + pOriginal * pOriginal) * MeV);
3020 vec[i]->SetMomentum(vec[i]->GetMomentum() * (pOriginal / testMomentum));
3021 }
3022 }
3023 }
3024
3025 void FullModelReactionDynamics::ProduceStrangeParticlePairs(G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> &vec,
3026 G4int &vecLen,
3027 const G4ReactionProduct &modifiedOriginal,
3028 const G4DynamicParticle *originalTarget,
3029 G4ReactionProduct ¤tParticle,
3030 G4ReactionProduct &targetParticle,
3031 G4bool &incidentHasChanged,
3032 G4bool &targetHasChanged) {
3033
3034
3035
3036
3037
3038
3039
3040
3041
3042 if (vecLen == 0)
3043 return;
3044
3045
3046
3047 if (currentParticle.GetMass() == 0.0 || targetParticle.GetMass() == 0.0)
3048 return;
3049
3050 const G4double etOriginal = modifiedOriginal.GetTotalEnergy() / GeV;
3051 const G4double mOriginal = modifiedOriginal.GetDefinition()->GetPDGMass() / GeV;
3052 G4double targetMass = originalTarget->GetDefinition()->GetPDGMass() / GeV;
3053 G4double centerofmassEnergy =
3054 std::sqrt(mOriginal * mOriginal + targetMass * targetMass + 2.0 * targetMass * etOriginal);
3055 G4double currentMass = currentParticle.GetMass() / GeV;
3056 G4double availableEnergy = centerofmassEnergy - (targetMass + currentMass);
3057 if (availableEnergy <= 1.0)
3058 return;
3059
3060 G4ParticleDefinition *aProton = G4Proton::Proton();
3061 G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
3062 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
3063 G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
3064 G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
3065 G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
3066 G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
3067 G4ParticleDefinition *anAntiSigmaMinus = G4AntiSigmaMinus::AntiSigmaMinus();
3068 G4ParticleDefinition *anAntiSigmaPlus = G4AntiSigmaPlus::AntiSigmaPlus();
3069 G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
3070 G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
3071 G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
3072 G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
3073 G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
3074 G4ParticleDefinition *aLambda = G4Lambda::Lambda();
3075 G4ParticleDefinition *anAntiLambda = G4AntiLambda::AntiLambda();
3076
3077 const G4double protonMass = aProton->GetPDGMass() / GeV;
3078 const G4double sigmaMinusMass = aSigmaMinus->GetPDGMass() / GeV;
3079
3080
3081
3082 const G4double avrs[] = {3., 4., 5., 6., 7., 8., 9., 10., 20., 30., 40., 50.};
3083
3084 G4int ibin, i3, i4;
3085 G4double avk, avy, avn, ran;
3086 G4int i = 1;
3087 while ((i < 12) && (centerofmassEnergy > avrs[i]))
3088 ++i;
3089 if (i == 12)
3090 ibin = 11;
3091 else
3092 ibin = i;
3093
3094
3095
3096
3097 if (vecLen == 1)
3098 {
3099 i3 = 0;
3100 i4 = 1;
3101 } else
3102 {
3103 G4double ran = G4UniformRand();
3104 while (ran == 1.0)
3105 ran = G4UniformRand();
3106 i4 = i3 = G4int(vecLen * ran);
3107 while (i3 == i4) {
3108 ran = G4UniformRand();
3109 while (ran == 1.0)
3110 ran = G4UniformRand();
3111 i4 = G4int(vecLen * ran);
3112 }
3113 }
3114
3115
3116
3117 const G4double avkkb[] = {0.0015, 0.005, 0.012, 0.0285, 0.0525, 0.075, 0.0975, 0.123, 0.28, 0.398, 0.495, 0.573};
3118 const G4double avky[] = {0.005, 0.03, 0.064, 0.095, 0.115, 0.13, 0.145, 0.155, 0.20, 0.205, 0.210, 0.212};
3119 const G4double avnnb[] = {0.00001, 0.0001, 0.0006, 0.0025, 0.01, 0.02, 0.04, 0.05, 0.12, 0.15, 0.18, 0.20};
3120
3121 avk = (std::log(avkkb[ibin]) - std::log(avkkb[ibin - 1])) * (centerofmassEnergy - avrs[ibin - 1]) /
3122 (avrs[ibin] - avrs[ibin - 1]) +
3123 std::log(avkkb[ibin - 1]);
3124 avk = std::exp(avk);
3125
3126 avy = (std::log(avky[ibin]) - std::log(avky[ibin - 1])) * (centerofmassEnergy - avrs[ibin - 1]) /
3127 (avrs[ibin] - avrs[ibin - 1]) +
3128 std::log(avky[ibin - 1]);
3129 avy = std::exp(avy);
3130
3131 avn = (std::log(avnnb[ibin]) - std::log(avnnb[ibin - 1])) * (centerofmassEnergy - avrs[ibin - 1]) /
3132 (avrs[ibin] - avrs[ibin - 1]) +
3133 std::log(avnnb[ibin - 1]);
3134 avn = std::exp(avn);
3135
3136 if (avk + avy + avn <= 0.0)
3137 return;
3138
3139 if (currentMass < protonMass)
3140 avy /= 2.0;
3141 if (targetMass < protonMass)
3142 avy = 0.0;
3143 avy += avk + avn;
3144 avk += avn;
3145 ran = G4UniformRand();
3146 if (ran < avn) {
3147 if (availableEnergy < 2.0)
3148 return;
3149 if (vecLen == 1)
3150 {
3151 G4ReactionProduct *p1 = new G4ReactionProduct;
3152 if (G4UniformRand() < 0.5) {
3153 vec[0]->SetDefinition(aNeutron);
3154 p1->SetDefinition(anAntiNeutron);
3155 (G4UniformRand() < 0.5) ? p1->SetSide(-1) : p1->SetSide(1);
3156 vec[0]->SetMayBeKilled(false);
3157 p1->SetMayBeKilled(false);
3158 } else {
3159 vec[0]->SetDefinition(aProton);
3160 p1->SetDefinition(anAntiProton);
3161 (G4UniformRand() < 0.5) ? p1->SetSide(-1) : p1->SetSide(1);
3162 vec[0]->SetMayBeKilled(false);
3163 p1->SetMayBeKilled(false);
3164 }
3165 vec.SetElement(vecLen++, p1);
3166
3167 } else {
3168 if (G4UniformRand() < 0.5) {
3169 vec[i3]->SetDefinition(aNeutron);
3170 vec[i4]->SetDefinition(anAntiNeutron);
3171 vec[i3]->SetMayBeKilled(false);
3172 vec[i4]->SetMayBeKilled(false);
3173 } else {
3174 vec[i3]->SetDefinition(aProton);
3175 vec[i4]->SetDefinition(anAntiProton);
3176 vec[i3]->SetMayBeKilled(false);
3177 vec[i4]->SetMayBeKilled(false);
3178 }
3179 }
3180 } else if (ran < avk) {
3181 if (availableEnergy < 1.0)
3182 return;
3183
3184 const G4double kkb[] = {0.2500, 0.3750, 0.5000, 0.5625, 0.6250, 0.6875, 0.7500, 0.8750, 1.000};
3185 const G4int ipakkb1[] = {10, 10, 10, 11, 11, 12, 12, 11, 12};
3186 const G4int ipakkb2[] = {13, 11, 12, 11, 12, 11, 12, 13, 13};
3187 ran = G4UniformRand();
3188 i = 0;
3189 while ((i < 9) && (ran >= kkb[i]))
3190 ++i;
3191 if (i == 9)
3192 return;
3193
3194
3195
3196
3197 switch (ipakkb1[i]) {
3198 case 10:
3199 vec[i3]->SetDefinition(aKaonPlus);
3200 vec[i3]->SetMayBeKilled(false);
3201 break;
3202 case 11:
3203 vec[i3]->SetDefinition(aKaonZS);
3204 vec[i3]->SetMayBeKilled(false);
3205 break;
3206 case 12:
3207 vec[i3]->SetDefinition(aKaonZL);
3208 vec[i3]->SetMayBeKilled(false);
3209 break;
3210 }
3211 if (vecLen == 1)
3212 {
3213 G4ReactionProduct *p1 = new G4ReactionProduct;
3214 switch (ipakkb2[i]) {
3215 case 11:
3216 p1->SetDefinition(aKaonZS);
3217 p1->SetMayBeKilled(false);
3218 break;
3219 case 12:
3220 p1->SetDefinition(aKaonZL);
3221 p1->SetMayBeKilled(false);
3222 break;
3223 case 13:
3224 p1->SetDefinition(aKaonMinus);
3225 p1->SetMayBeKilled(false);
3226 break;
3227 }
3228 (G4UniformRand() < 0.5) ? p1->SetSide(-1) : p1->SetSide(1);
3229 vec.SetElement(vecLen++, p1);
3230
3231 } else
3232 {
3233 switch (ipakkb2[i]) {
3234 case 11:
3235 vec[i4]->SetDefinition(aKaonZS);
3236 vec[i4]->SetMayBeKilled(false);
3237 break;
3238 case 12:
3239 vec[i4]->SetDefinition(aKaonZL);
3240 vec[i4]->SetMayBeKilled(false);
3241 break;
3242 case 13:
3243 vec[i4]->SetDefinition(aKaonMinus);
3244 vec[i4]->SetMayBeKilled(false);
3245 break;
3246 }
3247 }
3248 } else if (ran < avy) {
3249 if (availableEnergy < 1.6)
3250 return;
3251
3252 const G4double ky[] = {0.200, 0.300, 0.400, 0.550, 0.625, 0.700, 0.800, 0.850, 0.900, 0.950, 0.975, 1.000};
3253 const G4int ipaky1[] = {18, 18, 18, 20, 20, 20, 21, 21, 21, 22, 22, 22};
3254 const G4int ipaky2[] = {10, 11, 12, 10, 11, 12, 10, 11, 12, 10, 11, 12};
3255 const G4int ipakyb1[] = {19, 19, 19, 23, 23, 23, 24, 24, 24, 25, 25, 25};
3256 const G4int ipakyb2[] = {13, 12, 11, 13, 12, 11, 13, 12, 11, 13, 12, 11};
3257 ran = G4UniformRand();
3258 i = 0;
3259 while ((i < 12) && (ran > ky[i]))
3260 ++i;
3261 if (i == 12)
3262 return;
3263 if ((currentMass < protonMass) || (G4UniformRand() < 0.5)) {
3264
3265
3266
3267
3268
3269 switch (ipaky1[i]) {
3270 case 18:
3271 targetParticle.SetDefinition(aLambda);
3272 break;
3273 case 20:
3274 targetParticle.SetDefinition(aSigmaPlus);
3275 break;
3276 case 21:
3277 targetParticle.SetDefinition(aSigmaZero);
3278 break;
3279 case 22:
3280 targetParticle.SetDefinition(aSigmaMinus);
3281 break;
3282 }
3283 targetHasChanged = true;
3284 switch (ipaky2[i]) {
3285 case 10:
3286 vec[i3]->SetDefinition(aKaonPlus);
3287 vec[i3]->SetMayBeKilled(false);
3288 break;
3289 case 11:
3290 vec[i3]->SetDefinition(aKaonZS);
3291 vec[i3]->SetMayBeKilled(false);
3292 break;
3293 case 12:
3294 vec[i3]->SetDefinition(aKaonZL);
3295 vec[i3]->SetMayBeKilled(false);
3296 break;
3297 }
3298 } else
3299 {
3300 if ((currentParticle.GetDefinition() == anAntiProton) || (currentParticle.GetDefinition() == anAntiNeutron) ||
3301 (currentParticle.GetDefinition() == anAntiLambda) || (currentMass > sigmaMinusMass)) {
3302 switch (ipakyb1[i]) {
3303 case 19:
3304 currentParticle.SetDefinitionAndUpdateE(anAntiLambda);
3305 break;
3306 case 23:
3307 currentParticle.SetDefinitionAndUpdateE(anAntiSigmaPlus);
3308 break;
3309 case 24:
3310 currentParticle.SetDefinitionAndUpdateE(anAntiSigmaZero);
3311 break;
3312 case 25:
3313 currentParticle.SetDefinitionAndUpdateE(anAntiSigmaMinus);
3314 break;
3315 }
3316 incidentHasChanged = true;
3317 switch (ipakyb2[i]) {
3318 case 11:
3319 vec[i3]->SetDefinition(aKaonZS);
3320 vec[i3]->SetMayBeKilled(false);
3321 break;
3322 case 12:
3323 vec[i3]->SetDefinition(aKaonZL);
3324 vec[i3]->SetMayBeKilled(false);
3325 break;
3326 case 13:
3327 vec[i3]->SetDefinition(aKaonMinus);
3328 vec[i3]->SetMayBeKilled(false);
3329 break;
3330 }
3331 } else {
3332 switch (ipaky1[i]) {
3333 case 18:
3334 currentParticle.SetDefinitionAndUpdateE(aLambda);
3335 break;
3336 case 20:
3337 currentParticle.SetDefinitionAndUpdateE(aSigmaPlus);
3338 break;
3339 case 21:
3340 currentParticle.SetDefinitionAndUpdateE(aSigmaZero);
3341 break;
3342 case 22:
3343 currentParticle.SetDefinitionAndUpdateE(aSigmaMinus);
3344 break;
3345 }
3346 incidentHasChanged = true;
3347 switch (ipaky2[i]) {
3348 case 10:
3349 vec[i3]->SetDefinition(aKaonPlus);
3350 vec[i3]->SetMayBeKilled(false);
3351 break;
3352 case 11:
3353 vec[i3]->SetDefinition(aKaonZS);
3354 vec[i3]->SetMayBeKilled(false);
3355 break;
3356 case 12:
3357 vec[i3]->SetDefinition(aKaonZL);
3358 vec[i3]->SetMayBeKilled(false);
3359 break;
3360 }
3361 }
3362 }
3363 } else
3364 return;
3365
3366
3367
3368
3369
3370
3371
3372
3373
3374 currentMass = currentParticle.GetMass() / GeV;
3375 targetMass = targetParticle.GetMass() / GeV;
3376
3377 G4double energyCheck = centerofmassEnergy - (currentMass + targetMass);
3378 for (i = 0; i < vecLen; ++i) {
3379 energyCheck -= vec[i]->GetMass() / GeV;
3380 if (energyCheck < 0.0)
3381 {
3382 vecLen = std::max(0, --i);
3383 G4int j;
3384 for (j = i; j < vecLen; j++)
3385 delete vec[j];
3386 break;
3387 }
3388 }
3389 return;
3390 }
3391
3392 void FullModelReactionDynamics::NuclearReaction(G4FastVector<G4ReactionProduct, 4> &vec,
3393 G4int &vecLen,
3394 const G4HadProjectile *originalIncident,
3395 const G4Nucleus &targetNucleus,
3396 const G4double theAtomicMass,
3397 const G4double *mass) {
3398
3399
3400
3401
3402 G4ParticleDefinition *aGamma = G4Gamma::Gamma();
3403 G4ParticleDefinition *aProton = G4Proton::Proton();
3404 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
3405 G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
3406 G4ParticleDefinition *aTriton = G4Triton::Triton();
3407 G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
3408
3409 const G4double aProtonMass = aProton->GetPDGMass() / MeV;
3410 const G4double aNeutronMass = aNeutron->GetPDGMass() / MeV;
3411 const G4double aDeuteronMass = aDeuteron->GetPDGMass() / MeV;
3412 const G4double aTritonMass = aTriton->GetPDGMass() / MeV;
3413 const G4double anAlphaMass = anAlpha->GetPDGMass() / MeV;
3414
3415 G4ReactionProduct currentParticle;
3416 currentParticle = *originalIncident;
3417
3418
3419
3420
3421
3422 G4double p = currentParticle.GetTotalMomentum();
3423 G4double pp = currentParticle.GetMomentum().mag();
3424 if (pp <= 0.001 * MeV) {
3425 G4double phinve = twopi * G4UniformRand();
3426 G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3427 currentParticle.SetMomentum(
3428 p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3429 } else
3430 currentParticle.SetMomentum(currentParticle.GetMomentum() * (p / pp));
3431
3432
3433
3434 G4double currentKinetic = currentParticle.GetKineticEnergy() / MeV;
3435 G4double currentMass = currentParticle.GetDefinition()->GetPDGMass() / MeV;
3436 G4double qv = currentKinetic + theAtomicMass + currentMass;
3437
3438 G4double qval[9];
3439 qval[0] = qv - mass[0];
3440 qval[1] = qv - mass[1] - aNeutronMass;
3441 qval[2] = qv - mass[2] - aProtonMass;
3442 qval[3] = qv - mass[3] - aDeuteronMass;
3443 qval[4] = qv - mass[4] - aTritonMass;
3444 qval[5] = qv - mass[5] - anAlphaMass;
3445 qval[6] = qv - mass[6] - aNeutronMass - aNeutronMass;
3446 qval[7] = qv - mass[7] - aNeutronMass - aProtonMass;
3447 qval[8] = qv - mass[8] - aProtonMass - aProtonMass;
3448
3449 if (currentParticle.GetDefinition() == aNeutron) {
3450 const G4double A = targetNucleus.GetN_asInt();
3451 if (G4UniformRand() > ((A - 1.0) / 230.0) * ((A - 1.0) / 230.0))
3452 qval[0] = 0.0;
3453 if (G4UniformRand() >= currentKinetic / 7.9254 * A)
3454 qval[2] = qval[3] = qval[4] = qval[5] = qval[8] = 0.0;
3455 } else
3456 qval[0] = 0.0;
3457
3458 G4int i;
3459 qv = 0.0;
3460 for (i = 0; i < 9; ++i) {
3461 if (mass[i] < 500.0 * MeV)
3462 qval[i] = 0.0;
3463 if (qval[i] < 0.0)
3464 qval[i] = 0.0;
3465 qv += qval[i];
3466 }
3467 G4double qv1 = 0.0;
3468 G4double ran = G4UniformRand();
3469 G4int index;
3470 for (index = 0; index < 9; ++index) {
3471 if (qval[index] > 0.0) {
3472 qv1 += qval[index] / qv;
3473 if (ran <= qv1)
3474 break;
3475 }
3476 }
3477 if (index == 9)
3478 {
3479 throw G4HadronicException(
3480 __FILE__,
3481 __LINE__,
3482 "FullModelReactionDynamics::NuclearReaction: inelastic reaction kinematically not possible");
3483 }
3484 G4double ke = currentParticle.GetKineticEnergy() / GeV;
3485 G4int nt = 2;
3486 if ((index >= 6) || (G4UniformRand() < std::min(0.5, ke * 10.0)))
3487 nt = 3;
3488
3489 G4ReactionProduct **v = new G4ReactionProduct *[3];
3490 v[0] = new G4ReactionProduct;
3491 v[1] = new G4ReactionProduct;
3492 v[2] = new G4ReactionProduct;
3493
3494 v[0]->SetMass(mass[index] * MeV);
3495 switch (index) {
3496 case 0:
3497 v[1]->SetDefinition(aGamma);
3498 v[2]->SetDefinition(aGamma);
3499 break;
3500 case 1:
3501 v[1]->SetDefinition(aNeutron);
3502 v[2]->SetDefinition(aGamma);
3503 break;
3504 case 2:
3505 v[1]->SetDefinition(aProton);
3506 v[2]->SetDefinition(aGamma);
3507 break;
3508 case 3:
3509 v[1]->SetDefinition(aDeuteron);
3510 v[2]->SetDefinition(aGamma);
3511 break;
3512 case 4:
3513 v[1]->SetDefinition(aTriton);
3514 v[2]->SetDefinition(aGamma);
3515 break;
3516 case 5:
3517 v[1]->SetDefinition(anAlpha);
3518 v[2]->SetDefinition(aGamma);
3519 break;
3520 case 6:
3521 v[1]->SetDefinition(aNeutron);
3522 v[2]->SetDefinition(aNeutron);
3523 break;
3524 case 7:
3525 v[1]->SetDefinition(aNeutron);
3526 v[2]->SetDefinition(aProton);
3527 break;
3528 case 8:
3529 v[1]->SetDefinition(aProton);
3530 v[2]->SetDefinition(aProton);
3531 break;
3532 }
3533
3534
3535
3536 G4ReactionProduct pseudo1;
3537 pseudo1.SetMass(theAtomicMass * MeV);
3538 pseudo1.SetTotalEnergy(theAtomicMass * MeV);
3539 G4ReactionProduct pseudo2 = currentParticle + pseudo1;
3540 pseudo2.SetMomentum(pseudo2.GetMomentum() * (-1.0));
3541
3542
3543
3544 G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> tempV;
3545 tempV.Initialize(nt);
3546 G4int tempLen = 0;
3547 tempV.SetElement(tempLen++, v[0]);
3548 tempV.SetElement(tempLen++, v[1]);
3549 if (nt == 3)
3550 tempV.SetElement(tempLen++, v[2]);
3551 G4bool constantCrossSection = true;
3552 GenerateNBodyEvent(pseudo2.GetMass() / MeV, constantCrossSection, tempV, tempLen);
3553 v[0]->Lorentz(*v[0], pseudo2);
3554 v[1]->Lorentz(*v[1], pseudo2);
3555 if (nt == 3)
3556 v[2]->Lorentz(*v[2], pseudo2);
3557
3558 G4bool particleIsDefined = false;
3559 if (v[0]->GetMass() / MeV - aProtonMass < 0.1) {
3560 v[0]->SetDefinition(aProton);
3561 particleIsDefined = true;
3562 } else if (v[0]->GetMass() / MeV - aNeutronMass < 0.1) {
3563 v[0]->SetDefinition(aNeutron);
3564 particleIsDefined = true;
3565 } else if (v[0]->GetMass() / MeV - aDeuteronMass < 0.1) {
3566 v[0]->SetDefinition(aDeuteron);
3567 particleIsDefined = true;
3568 } else if (v[0]->GetMass() / MeV - aTritonMass < 0.1) {
3569 v[0]->SetDefinition(aTriton);
3570 particleIsDefined = true;
3571 } else if (v[0]->GetMass() / MeV - anAlphaMass < 0.1) {
3572 v[0]->SetDefinition(anAlpha);
3573 particleIsDefined = true;
3574 }
3575 currentParticle.SetKineticEnergy(std::max(0.001, currentParticle.GetKineticEnergy() / MeV));
3576 p = currentParticle.GetTotalMomentum();
3577 pp = currentParticle.GetMomentum().mag();
3578 if (pp <= 0.001 * MeV) {
3579 G4double phinve = twopi * G4UniformRand();
3580 G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3581 currentParticle.SetMomentum(
3582 p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3583 } else
3584 currentParticle.SetMomentum(currentParticle.GetMomentum() * (p / pp));
3585
3586 if (particleIsDefined) {
3587 v[0]->SetKineticEnergy(std::max(0.001, 0.5 * G4UniformRand() * v[0]->GetKineticEnergy() / MeV));
3588 p = v[0]->GetTotalMomentum();
3589 pp = v[0]->GetMomentum().mag();
3590 if (pp <= 0.001 * MeV) {
3591 G4double phinve = twopi * G4UniformRand();
3592 G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3593 v[0]->SetMomentum(
3594 p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3595 } else
3596 v[0]->SetMomentum(v[0]->GetMomentum() * (p / pp));
3597 }
3598 if ((v[1]->GetDefinition() == aDeuteron) || (v[1]->GetDefinition() == aTriton) || (v[1]->GetDefinition() == anAlpha))
3599 v[1]->SetKineticEnergy(std::max(0.001, 0.5 * G4UniformRand() * v[1]->GetKineticEnergy() / MeV));
3600 else
3601 v[1]->SetKineticEnergy(std::max(0.001, v[1]->GetKineticEnergy() / MeV));
3602
3603 p = v[1]->GetTotalMomentum();
3604 pp = v[1]->GetMomentum().mag();
3605 if (pp <= 0.001 * MeV) {
3606 G4double phinve = twopi * G4UniformRand();
3607 G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3608 v[1]->SetMomentum(
3609 p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3610 } else
3611 v[1]->SetMomentum(v[1]->GetMomentum() * (p / pp));
3612
3613 if (nt == 3) {
3614 if ((v[2]->GetDefinition() == aDeuteron) || (v[2]->GetDefinition() == aTriton) ||
3615 (v[2]->GetDefinition() == anAlpha))
3616 v[2]->SetKineticEnergy(std::max(0.001, 0.5 * G4UniformRand() * v[2]->GetKineticEnergy() / MeV));
3617 else
3618 v[2]->SetKineticEnergy(std::max(0.001, v[2]->GetKineticEnergy() / MeV));
3619
3620 p = v[2]->GetTotalMomentum();
3621 pp = v[2]->GetMomentum().mag();
3622 if (pp <= 0.001 * MeV) {
3623 G4double phinve = twopi * G4UniformRand();
3624 G4double rthnve = std::acos(std::max(-1.0, std::min(1.0, -1.0 + 2.0 * G4UniformRand())));
3625 v[2]->SetMomentum(
3626 p * std::sin(rthnve) * std::cos(phinve), p * std::sin(rthnve) * std::sin(phinve), p * std::cos(rthnve));
3627 } else
3628 v[2]->SetMomentum(v[2]->GetMomentum() * (p / pp));
3629 }
3630 G4int del;
3631 for (del = 0; del < vecLen; del++)
3632 delete vec[del];
3633 vecLen = 0;
3634 if (particleIsDefined) {
3635 vec.SetElement(vecLen++, v[0]);
3636 } else {
3637 delete v[0];
3638 }
3639 vec.SetElement(vecLen++, v[1]);
3640 if (nt == 3) {
3641 vec.SetElement(vecLen++, v[2]);
3642 } else {
3643 delete v[2];
3644 }
3645 delete[] v;
3646 return;
3647 }
3648
3649