File indexing completed on 2024-11-22 02:45:26
0001 #include "G4ProcessManager.hh"
0002 #include "G4ParticleTable.hh"
0003 #include "G4HadronicException.hh"
0004
0005 #include "SimG4Core/CustomPhysics/interface/FullModelHadronicProcess.h"
0006 #include "SimG4Core/CustomPhysics/interface/G4ProcessHelper.h"
0007 #include "SimG4Core/CustomPhysics/interface/Decay3Body.h"
0008 #include "SimG4Core/CustomPhysics/interface/CustomPDGParser.h"
0009 #include "SimG4Core/CustomPhysics/interface/CustomParticle.h"
0010
0011 using namespace CLHEP;
0012
0013 FullModelHadronicProcess::FullModelHadronicProcess(G4ProcessHelper* aHelper, const G4String& processName)
0014 : G4VDiscreteProcess(processName), theHelper(aHelper) {}
0015
0016 FullModelHadronicProcess::~FullModelHadronicProcess() {}
0017
0018 G4bool FullModelHadronicProcess::IsApplicable(const G4ParticleDefinition& aP) {
0019 return theHelper->ApplicabilityTester(aP);
0020 }
0021
0022 G4double FullModelHadronicProcess::GetMicroscopicCrossSection(const G4DynamicParticle* aParticle,
0023 const G4Element* anElement,
0024 G4double aTemp) {
0025
0026 G4double InclXsec = theHelper->GetInclusiveCrossSection(aParticle, anElement);
0027
0028 return InclXsec;
0029 }
0030
0031 G4double FullModelHadronicProcess::GetMeanFreePath(const G4Track& aTrack, G4double, G4ForceCondition*) {
0032 G4Material* aMaterial = aTrack.GetMaterial();
0033 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
0034 G4double sigma = 0.0;
0035
0036 G4int nElements = aMaterial->GetNumberOfElements();
0037
0038 const G4double* theAtomicNumDensityVector = aMaterial->GetAtomicNumDensityVector();
0039 G4double aTemp = aMaterial->GetTemperature();
0040
0041 for (G4int i = 0; i < nElements; ++i) {
0042 G4double xSection = GetMicroscopicCrossSection(aParticle, (*aMaterial->GetElementVector())[i], aTemp);
0043 sigma += theAtomicNumDensityVector[i] * xSection;
0044 }
0045 G4double res = DBL_MAX;
0046 if (sigma > 0.0) {
0047 res = 1. / sigma;
0048 }
0049 return res;
0050 }
0051
0052 G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) {
0053
0054 const G4TouchableHandle& thisTouchable(aTrack.GetTouchableHandle());
0055
0056
0057 aParticleChange.Initialize(aTrack);
0058 const G4DynamicParticle* IncidentRhadron = aTrack.GetDynamicParticle();
0059 CustomParticle* CustomIncident = static_cast<CustomParticle*>(IncidentRhadron->GetDefinition());
0060 const G4ThreeVector& aPosition = aTrack.GetPosition();
0061
0062 const G4int theIncidentPDG = IncidentRhadron->GetDefinition()->GetPDGEncoding();
0063 G4ParticleTable* theParticleTable = G4ParticleTable::GetParticleTable();
0064 std::vector<G4ParticleDefinition*> theParticleDefinitions;
0065 G4bool IncidentSurvives = false;
0066 G4bool TargetSurvives = false;
0067 G4Nucleus targetNucleus(aTrack.GetMaterial());
0068 G4ParticleDefinition* outgoingRhadron = nullptr;
0069 G4ParticleDefinition* outgoingCloud = nullptr;
0070 G4ParticleDefinition* outgoingTarget = nullptr;
0071
0072 G4ThreeVector p_0 = IncidentRhadron->GetMomentum();
0073 G4double e_kin_0 = IncidentRhadron->GetKineticEnergy();
0074
0075
0076 G4DynamicParticle* cloudParticle = new G4DynamicParticle();
0077
0078
0079
0080
0081
0082
0083 cloudParticle->SetDefinition(CustomIncident->GetCloud());
0084
0085 if (cloudParticle->GetDefinition() == nullptr) {
0086 G4cout << "FullModelHadronicProcess::PostStepDoIt Definition of particle cloud not available!!" << G4endl;
0087 }
0088
0089
0090
0091
0092
0093 double scale = cloudParticle->GetDefinition()->GetPDGMass() / IncidentRhadron->GetDefinition()->GetPDGMass();
0094
0095 G4LorentzVector cloudMomentum(IncidentRhadron->GetMomentum() * scale, cloudParticle->GetDefinition()->GetPDGMass());
0096 G4LorentzVector gluinoMomentum(IncidentRhadron->GetMomentum() * (1. - scale),
0097 CustomIncident->GetSpectator()->GetPDGMass());
0098
0099
0100 G4LorentzVector FullRhadron4Momentum = IncidentRhadron->Get4Momentum();
0101 const G4LorentzVector& Cloud4Momentum = cloudMomentum;
0102
0103 cloudParticle->Set4Momentum(cloudMomentum);
0104
0105 G4DynamicParticle* OrgPart = cloudParticle;
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117 double E_0 = IncidentRhadron->GetKineticEnergy();
0118 G4double ek = OrgPart->GetKineticEnergy();
0119 G4double amas = OrgPart->GetDefinition()->GetPDGMass();
0120 G4ThreeVector dir = (OrgPart->GetMomentum()).unit();
0121 G4double tkin = targetNucleus.Cinema(ek);
0122 ek += tkin;
0123
0124
0125 tkin = targetNucleus.EvaporationEffects(ek);
0126 ek -= tkin;
0127
0128 if (ek + gluinoMomentum.e() - gluinoMomentum.m() <= 0.1 * MeV || ek <= 0.) {
0129
0130 G4cout << "Kinetic energy is sick" << G4endl;
0131 G4cout << "Full R-hadron: " << (ek + gluinoMomentum.e() - gluinoMomentum.m()) / MeV << " MeV" << G4endl;
0132 G4cout << "Quark system: " << ek / MeV << " MeV" << G4endl;
0133 aParticleChange.ProposeTrackStatus(fStopButAlive);
0134 return &aParticleChange;
0135 }
0136 OrgPart->SetKineticEnergy(ek);
0137 G4double p = std::sqrt(ek * (ek + 2 * amas));
0138 OrgPart->SetMomentum(dir * p);
0139
0140
0141 G4ParticleDefinition* aTarget;
0142 ReactionProduct rp = theHelper->GetFinalState(aTrack, aTarget);
0143 G4bool force2to2 = false;
0144
0145 while (rp.size() != 2 && force2to2) {
0146 rp = theHelper->GetFinalState(aTrack, aTarget);
0147 }
0148 G4int NumberOfSecondaries = rp.size();
0149
0150
0151
0152 G4LorentzVector Target4Momentum(0., 0., 0., aTarget->GetPDGMass());
0153
0154 G4LorentzVector psum_full = FullRhadron4Momentum + Target4Momentum;
0155 G4LorentzVector psum_cloud = Cloud4Momentum + Target4Momentum;
0156 G4ThreeVector trafo_full_cms = (-1) * psum_full.boostVector();
0157 G4ThreeVector trafo_cloud_cms = (-1) * psum_cloud.boostVector();
0158
0159
0160
0161 for (ReactionProduct::iterator it = rp.begin(); it != rp.end(); ++it) {
0162 G4ParticleDefinition* tempDef = theParticleTable->FindParticle(*it);
0163 CustomParticle* tempCust = dynamic_cast<CustomParticle*>(tempDef);
0164 if (tempDef == aTarget)
0165 TargetSurvives = true;
0166
0167
0168 if (tempCust != nullptr) {
0169 outgoingRhadron = tempDef;
0170
0171 outgoingCloud = tempCust->GetCloud();
0172 if (outgoingCloud == nullptr) {
0173 G4cout << "FullModelHadronicProcess::PostStepDoIt Definition of outgoing particle cloud not available!"
0174 << G4endl;
0175 }
0176
0177
0178
0179
0180 }
0181
0182 if (tempDef == G4Proton::Proton() || tempDef == G4Neutron::Neutron())
0183 outgoingTarget = tempDef;
0184 if (tempCust == nullptr && rp.size() == 2)
0185 outgoingTarget = tempDef;
0186 if (tempDef->GetPDGEncoding() == theIncidentPDG) {
0187 IncidentSurvives = true;
0188 } else {
0189 theParticleDefinitions.push_back(tempDef);
0190 }
0191 }
0192
0193 if (outgoingTarget == nullptr)
0194 outgoingTarget = theParticleTable->FindParticle(rp[1]);
0195
0196
0197
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207
0208 if (IncidentSurvives)
0209 NumberOfSecondaries--;
0210 aParticleChange.SetNumberOfSecondaries(NumberOfSecondaries);
0211
0212
0213
0214
0215
0216
0217 const G4HadProjectile* originalIncident = new G4HadProjectile(*OrgPart);
0218
0219
0220 G4HadProjectile* org2 = new G4HadProjectile(*OrgPart);
0221 G4LorentzRotation trans = org2->GetTrafoToLab();
0222 delete org2;
0223
0224
0225
0226 G4DynamicParticle* originalTarget = new G4DynamicParticle;
0227 originalTarget->SetDefinition(aTarget);
0228
0229 G4ReactionProduct targetParticle(aTarget);
0230
0231 G4ReactionProduct currentParticle(const_cast<G4ParticleDefinition*>(originalIncident->GetDefinition()));
0232 currentParticle.SetMomentum(originalIncident->Get4Momentum().vect());
0233 currentParticle.SetKineticEnergy(originalIncident->GetKineticEnergy());
0234
0235
0236
0237
0238
0239
0240
0241
0242 G4ReactionProduct modifiedOriginal = currentParticle;
0243
0244 currentParticle.SetSide(1);
0245 targetParticle.SetSide(-1);
0246 G4bool incidentHasChanged = false;
0247 if (!IncidentSurvives)
0248 incidentHasChanged = true;
0249 G4bool targetHasChanged = false;
0250 if (!TargetSurvives)
0251 targetHasChanged = true;
0252 G4bool quasiElastic = false;
0253 if (rp.size() == 2)
0254 quasiElastic = true;
0255 G4FastVector<G4ReactionProduct, MYGHADLISTSIZE> vec;
0256 G4int vecLen = 0;
0257 vec.Initialize(0);
0258
0259
0260
0261
0262 for (G4int i = 0; i != NumberOfSecondaries; i++) {
0263 if (theParticleDefinitions[i] != aTarget && theParticleDefinitions[i] != originalIncident->GetDefinition() &&
0264 theParticleDefinitions[i] != outgoingRhadron && theParticleDefinitions[i] != outgoingTarget) {
0265 G4ReactionProduct* pa = new G4ReactionProduct;
0266 pa->SetDefinition(theParticleDefinitions[i]);
0267 (G4UniformRand() < 0.5) ? pa->SetSide(-1) : pa->SetSide(1);
0268 vec.SetElement(vecLen++, pa);
0269 }
0270 }
0271
0272 if (incidentHasChanged)
0273 currentParticle.SetDefinitionAndUpdateE(outgoingCloud);
0274 if (incidentHasChanged)
0275 modifiedOriginal.SetDefinition(
0276 outgoingCloud);
0277 if (targetHasChanged)
0278 targetParticle.SetDefinitionAndUpdateE(outgoingTarget);
0279
0280
0281
0282
0283
0284
0285
0286
0287 CalculateMomenta(vec,
0288 vecLen,
0289 originalIncident,
0290 originalTarget,
0291 modifiedOriginal,
0292 targetNucleus,
0293 currentParticle,
0294 targetParticle,
0295 incidentHasChanged,
0296 targetHasChanged,
0297 quasiElastic);
0298
0299
0300
0301 G4String cPname = currentParticle.GetDefinition()->GetParticleName();
0302
0303
0304
0305
0306
0307
0308
0309
0310
0311
0312
0313
0314
0315
0316
0317
0318
0319
0320
0321
0322
0323
0324
0325
0326
0327
0328 aParticleChange.SetNumberOfSecondaries(vecLen + NumberOfSecondaries);
0329 G4double e_kin = 0;
0330 G4LorentzVector cloud_p4_new;
0331
0332
0333
0334
0335
0336
0337
0338
0339
0340
0341
0342 cloud_p4_new.setVectM(currentParticle.GetMomentum(), currentParticle.GetMass());
0343 cloud_p4_new *= trans;
0344
0345 G4LorentzVector cloud_p4_old_full = Cloud4Momentum;
0346 cloud_p4_old_full.boost(trafo_full_cms);
0347 G4LorentzVector cloud_p4_old_cloud = Cloud4Momentum;
0348 cloud_p4_old_cloud.boost(trafo_cloud_cms);
0349 G4LorentzVector cloud_p4_full = cloud_p4_new;
0350 cloud_p4_full.boost(trafo_full_cms);
0351 G4LorentzVector cloud_p4_cloud = cloud_p4_new;
0352 cloud_p4_cloud.boost(trafo_cloud_cms);
0353
0354 G4LorentzVector p_g_cms = gluinoMomentum;
0355 p_g_cms.boost(trafo_full_cms);
0356
0357 G4LorentzVector p4_new(cloud_p4_new.v() + gluinoMomentum.v(), outgoingRhadron->GetPDGMass());
0358
0359
0360
0361 G4ThreeVector p_new = p4_new.vect();
0362
0363 aParticleChange.ProposeLocalEnergyDeposit((p4_new - cloud_p4_new - gluinoMomentum).m());
0364
0365 if (incidentHasChanged) {
0366 G4DynamicParticle* p0 = new G4DynamicParticle;
0367 p0->SetDefinition(outgoingRhadron);
0368 p0->SetMomentum(p_new);
0369
0370
0371 G4Track* Track0 = new G4Track(p0, aTrack.GetGlobalTime(), aPosition);
0372 Track0->SetTouchableHandle(thisTouchable);
0373 aParticleChange.AddSecondary(Track0);
0374
0375
0376
0377
0378
0379 if (p0->GetKineticEnergy() > e_kin_0) {
0380 G4cout << "ALAAAAARM!!! (incident changed from " << IncidentRhadron->GetDefinition()->GetParticleName() << " to "
0381 << p0->GetDefinition()->GetParticleName() << ")" << G4endl;
0382 G4cout << "Energy loss: " << (e_kin_0 - p0->GetKineticEnergy()) / GeV << " GeV (should be positive)" << G4endl;
0383
0384 if (rp.size() != 3)
0385 G4cout << "DOUBLE ALAAAAARM!!!" << G4endl;
0386 }
0387
0388
0389 if (std::abs(p0->GetKineticEnergy() - e_kin_0) > 100 * GeV) {
0390 G4cout << "Diff. too big" << G4endl;
0391 }
0392 aParticleChange.ProposeTrackStatus(fStopAndKill);
0393 } else {
0394 G4double p = p_new.mag();
0395 if (p > DBL_MIN)
0396 aParticleChange.ProposeMomentumDirection(p_new.x() / p, p_new.y() / p, p_new.z() / p);
0397 else
0398 aParticleChange.ProposeMomentumDirection(1.0, 0.0, 0.0);
0399 }
0400
0401
0402 if (targetParticle.GetMass() > 0.0)
0403 {
0404 G4DynamicParticle* p1 = new G4DynamicParticle;
0405 p1->SetDefinition(targetParticle.GetDefinition());
0406
0407 G4ThreeVector momentum = targetParticle.GetMomentum();
0408 momentum = momentum.rotate(cache, what);
0409 p1->SetMomentum(momentum);
0410 p1->SetMomentum((trans * p1->Get4Momentum()).vect());
0411 G4Track* Track1 = new G4Track(p1, aTrack.GetGlobalTime(), aPosition);
0412 Track1->SetTouchableHandle(thisTouchable);
0413 aParticleChange.AddSecondary(Track1);
0414 }
0415 G4DynamicParticle* pa;
0416
0417
0418
0419
0420
0421 for (int i = 0; i < vecLen; ++i) {
0422 pa = new G4DynamicParticle();
0423 pa->SetDefinition(vec[i]->GetDefinition());
0424 pa->SetMomentum(vec[i]->GetMomentum());
0425 pa->Set4Momentum(trans * (pa->Get4Momentum()));
0426 G4ThreeVector pvec = pa->GetMomentum();
0427 G4Track* Trackn = new G4Track(pa, aTrack.GetGlobalTime(), aPosition);
0428 Trackn->SetTouchableHandle(thisTouchable);
0429 aParticleChange.AddSecondary(Trackn);
0430
0431 delete vec[i];
0432 }
0433
0434
0435 const G4DynamicParticle* theRhadron = FindRhadron(&aParticleChange);
0436
0437 if (theRhadron != nullptr || IncidentSurvives) {
0438 double E_new;
0439 if (IncidentSurvives) {
0440 E_new = e_kin;
0441 } else {
0442 E_new = theRhadron->GetKineticEnergy();
0443 if (CustomPDGParser::s_isRMeson(theRhadron->GetDefinition()->GetPDGEncoding()) !=
0444 CustomPDGParser::s_isRMeson(theIncidentPDG) ||
0445 CustomPDGParser::s_isMesonino(theRhadron->GetDefinition()->GetPDGEncoding()) !=
0446 CustomPDGParser::s_isMesonino(theIncidentPDG)) {
0447 G4cout << "Rm: " << CustomPDGParser::s_isRMeson(theRhadron->GetDefinition()->GetPDGEncoding())
0448 << " vs: " << CustomPDGParser::s_isRMeson(theIncidentPDG) << G4endl;
0449 G4cout << "Sm: " << CustomPDGParser::s_isMesonino(theRhadron->GetDefinition()->GetPDGEncoding())
0450 << " vs: " << CustomPDGParser::s_isMesonino(theIncidentPDG) << G4endl;
0451 }
0452 }
0453
0454
0455 G4LorentzVector p4_old_full = FullRhadron4Momentum;
0456 p4_old_full.boost(trafo_full_cms);
0457 G4LorentzVector p4_old_cloud = FullRhadron4Momentum;
0458 p4_old_cloud.boost(trafo_cloud_cms);
0459 G4LorentzVector p4_full = p4_new;
0460
0461 p4_full = p4_full.boost(trafo_full_cms);
0462
0463 G4LorentzVector p4_cloud = p4_new;
0464 p4_cloud.boost(trafo_cloud_cms);
0465
0466 G4double AbsDeltaE = E_0 - E_new;
0467
0468 if (AbsDeltaE > 10 * GeV) {
0469 G4cout << "Energy loss larger than 10 GeV..." << G4endl;
0470 G4cout << "E_0: " << E_0 / GeV << " GeV" << G4endl;
0471 G4cout << "E_new: " << E_new / GeV << " GeV" << G4endl;
0472 G4cout << "Gamma: " << IncidentRhadron->GetTotalEnergy() / IncidentRhadron->GetDefinition()->GetPDGMass()
0473 << G4endl;
0474 G4cout << "x: " << aPosition.x() / cm << " cm" << G4endl;
0475 }
0476 }
0477 delete originalIncident;
0478 delete originalTarget;
0479
0480
0481
0482
0483 ClearNumberOfInteractionLengthLeft();
0484
0485 return &aParticleChange;
0486 }
0487
0488 void FullModelHadronicProcess::CalculateMomenta(
0489 G4FastVector<G4ReactionProduct, MYGHADLISTSIZE>& vec,
0490 G4int& vecLen,
0491 const G4HadProjectile* originalIncident,
0492 const G4DynamicParticle* originalTarget,
0493 G4ReactionProduct& modifiedOriginal,
0494 G4Nucleus& targetNucleus,
0495 G4ReactionProduct& currentParticle,
0496 G4ReactionProduct& targetParticle,
0497 G4bool& incidentHasChanged,
0498 G4bool& targetHasChanged,
0499 G4bool quasiElastic) {
0500 FullModelReactionDynamics theReactionDynamics;
0501
0502 cache = 0;
0503 what = originalIncident->Get4Momentum().vect();
0504
0505 if (quasiElastic) {
0506
0507 theReactionDynamics.TwoBody(
0508 vec, vecLen, modifiedOriginal, originalTarget, currentParticle, targetParticle, targetNucleus, targetHasChanged);
0509
0510 return;
0511 }
0512
0513
0514 G4ReactionProduct leadingStrangeParticle;
0515 G4bool leadFlag = MarkLeadingStrangeParticle(currentParticle, targetParticle, leadingStrangeParticle);
0516
0517
0518
0519
0520 G4bool finishedGenXPt = false;
0521 G4bool annihilation = false;
0522 if (originalIncident->GetDefinition()->GetPDGEncoding() < 0 && currentParticle.GetMass() == 0.0 &&
0523 targetParticle.GetMass() == 0.0) {
0524
0525 annihilation = true;
0526 G4double ekcor = 1.0;
0527 G4double ek = originalIncident->GetKineticEnergy();
0528 G4double ekOrg = ek;
0529
0530 const G4double tarmas = originalTarget->GetDefinition()->GetPDGMass();
0531 if (ek > 1.0 * GeV)
0532 ekcor = 1. / (ek / GeV);
0533 const G4double atomicWeight = G4double(targetNucleus.GetN_asInt());
0534 ek = 2 * tarmas + ek * (1. + ekcor / atomicWeight);
0535
0536 G4double tkin = targetNucleus.Cinema(ek);
0537
0538 ekOrg += tkin;
0539 modifiedOriginal.SetKineticEnergy(ekOrg);
0540 }
0541
0542 const G4double twsup[] = {1.0, 0.7, 0.5, 0.3, 0.2, 0.1};
0543 G4double rand1 = G4UniformRand();
0544 G4double rand2 = G4UniformRand();
0545 if ((annihilation || (vecLen >= 6) || (modifiedOriginal.GetKineticEnergy() / GeV >= 1.0)) &&
0546 (((originalIncident->GetDefinition() == G4KaonPlus::KaonPlus()) ||
0547 (originalIncident->GetDefinition() == G4KaonMinus::KaonMinus()) ||
0548 (originalIncident->GetDefinition() == G4KaonZeroLong::KaonZeroLong()) ||
0549 (originalIncident->GetDefinition() == G4KaonZeroShort::KaonZeroShort())) &&
0550 ((rand1 < 0.5) || (rand2 > twsup[vecLen]))))
0551 finishedGenXPt = theReactionDynamics.GenerateXandPt(vec,
0552 vecLen,
0553 modifiedOriginal,
0554 originalIncident,
0555 currentParticle,
0556 targetParticle,
0557 targetNucleus,
0558 incidentHasChanged,
0559 targetHasChanged,
0560 leadFlag,
0561 leadingStrangeParticle);
0562 if (finishedGenXPt) {
0563 Rotate(vec, vecLen);
0564 return;
0565 }
0566
0567 G4bool finishedTwoClu = false;
0568 if (modifiedOriginal.GetTotalMomentum() / MeV < 1.0) {
0569 for (G4int i = 0; i < vecLen; i++)
0570 delete vec[i];
0571 vecLen = 0;
0572 } else {
0573 theReactionDynamics.SuppressChargedPions(vec,
0574 vecLen,
0575 modifiedOriginal,
0576 currentParticle,
0577 targetParticle,
0578 targetNucleus,
0579 incidentHasChanged,
0580 targetHasChanged);
0581
0582 try {
0583 finishedTwoClu = theReactionDynamics.TwoCluster(vec,
0584 vecLen,
0585 modifiedOriginal,
0586 originalIncident,
0587 currentParticle,
0588 targetParticle,
0589 targetNucleus,
0590 incidentHasChanged,
0591 targetHasChanged,
0592 leadFlag,
0593 leadingStrangeParticle);
0594 } catch (G4HadronicException& aR) {
0595 G4ExceptionDescription ed;
0596 aR.Report(ed);
0597 G4Exception("FullModelHadronicProcess::CalculateMomenta", "had066", FatalException, ed);
0598 }
0599 }
0600 if (finishedTwoClu) {
0601 Rotate(vec, vecLen);
0602 return;
0603 }
0604
0605
0606
0607
0608
0609
0610
0611
0612
0613
0614
0615
0616 theReactionDynamics.TwoBody(
0617 vec, vecLen, modifiedOriginal, originalTarget, currentParticle, targetParticle, targetNucleus, targetHasChanged);
0618 }
0619
0620 G4bool FullModelHadronicProcess::MarkLeadingStrangeParticle(const G4ReactionProduct& currentParticle,
0621 const G4ReactionProduct& targetParticle,
0622 G4ReactionProduct& leadParticle) {
0623
0624
0625
0626
0627
0628 G4bool lead = false;
0629 if ((currentParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
0630 (currentParticle.GetDefinition() != G4Proton::Proton()) &&
0631 (currentParticle.GetDefinition() != G4Neutron::Neutron())) {
0632 lead = true;
0633 leadParticle = currentParticle;
0634 } else if ((targetParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
0635 (targetParticle.GetDefinition() != G4Proton::Proton()) &&
0636 (targetParticle.GetDefinition() != G4Neutron::Neutron())) {
0637 lead = true;
0638 leadParticle = targetParticle;
0639 }
0640 return lead;
0641 }
0642
0643 void FullModelHadronicProcess::Rotate(G4FastVector<G4ReactionProduct, MYGHADLISTSIZE>& vec, G4int& vecLen) {
0644 G4double rotation = 2. * pi * G4UniformRand();
0645 cache = rotation;
0646 G4int i;
0647 for (i = 0; i < vecLen; ++i) {
0648 G4ThreeVector momentum = vec[i]->GetMomentum();
0649 momentum = momentum.rotate(rotation, what);
0650 vec[i]->SetMomentum(momentum);
0651 }
0652 }
0653
0654 const G4DynamicParticle* FullModelHadronicProcess::FindRhadron(G4ParticleChange* aParticleChange) {
0655 G4int nsec = aParticleChange->GetNumberOfSecondaries();
0656 if (nsec == 0)
0657 return nullptr;
0658 int i = 0;
0659 G4bool found = false;
0660 while (i != nsec && !found) {
0661
0662
0663 if (dynamic_cast<CustomParticle*>(aParticleChange->GetSecondary(i)->GetDynamicParticle()->GetDefinition()) !=
0664 nullptr)
0665 found = true;
0666 i++;
0667 }
0668 i--;
0669 if (found)
0670 return aParticleChange->GetSecondary(i)->GetDynamicParticle();
0671 return nullptr;
0672 }