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