File indexing completed on 2024-05-10 02:21:24
0001
0002 #include "SimG4Core/CustomPhysics/interface/CMSSIMPInelasticProcess.h"
0003 #include "SimG4Core/CustomPhysics/interface/CMSSIMPInelasticXS.h"
0004 #include "SimG4Core/CustomPhysics/interface/CMSSIMP.h"
0005
0006 #include "G4Types.hh"
0007 #include <CLHEP/Units/SystemOfUnits.h>
0008 #include "G4HadProjectile.hh"
0009 #include "G4ElementVector.hh"
0010 #include "G4Track.hh"
0011 #include "G4Step.hh"
0012 #include "G4Element.hh"
0013 #include "G4ParticleChange.hh"
0014 #include "G4NucleiProperties.hh"
0015 #include "G4Nucleus.hh"
0016
0017 #include "G4HadronicException.hh"
0018 #include "G4HadronicProcessStore.hh"
0019 #include "G4HadronicInteraction.hh"
0020
0021 #include "G4ParticleDefinition.hh"
0022
0023
0024 CMSSIMPInelasticProcess::CMSSIMPInelasticProcess(const G4String& processName)
0025 : G4HadronicProcess(processName, fHadronic) {
0026 AddDataSet(new CMSSIMPInelasticXS());
0027 theParticle = CMSSIMP::SIMP();
0028 }
0029
0030 CMSSIMPInelasticProcess::~CMSSIMPInelasticProcess() {}
0031
0032 G4bool CMSSIMPInelasticProcess::IsApplicable(const G4ParticleDefinition& aP) {
0033 return theParticle->GetParticleType() == aP.GetParticleType();
0034 }
0035
0036 G4VParticleChange* CMSSIMPInelasticProcess::PostStepDoIt(const G4Track& aTrack, const G4Step&) {
0037
0038 theTotalResult->Clear();
0039 theTotalResult->Initialize(aTrack);
0040 theTotalResult->ProposeWeight(aTrack.GetWeight());
0041 if (aTrack.GetTrackStatus() != fAlive) {
0042 return theTotalResult;
0043 }
0044
0045
0046
0047 G4DynamicParticle* aParticle = const_cast<G4DynamicParticle*>(aTrack.GetDynamicParticle());
0048
0049
0050 aParticle->SetPDGcode(2112);
0051 aParticle->SetDefinition(G4Neutron::Neutron());
0052
0053 G4Nucleus* target = GetTargetNucleusPointer();
0054 const G4Material* aMaterial = aTrack.GetMaterial();
0055 const G4Element* anElement = GetCrossSectionDataStore()->SampleZandA(aParticle, aMaterial, *target);
0056
0057
0058
0059 if (aTrack.GetTrackStatus() != fAlive && aTrack.GetTrackStatus() != fSuspend) {
0060 if (aTrack.GetTrackStatus() == fStopAndKill || aTrack.GetTrackStatus() == fKillTrackAndSecondaries ||
0061 aTrack.GetTrackStatus() == fPostponeToNextEvent) {
0062 G4ExceptionDescription ed;
0063 ed << "CMSSIMPInelasticProcess: track in unusable state - " << aTrack.GetTrackStatus() << G4endl;
0064 ed << "CMSSIMPInelasticProcess: returning unchanged track " << G4endl;
0065 DumpState(aTrack, "PostStepDoIt", ed);
0066 G4Exception("CMSSIMPInelasticProcess::PostStepDoIt", "had004", JustWarning, ed);
0067 }
0068
0069 return theTotalResult;
0070 }
0071
0072
0073 thePro.Initialise(aTrack);
0074 G4HadronicInteraction* anInteraction = GetHadronicInteractionList()[0];
0075
0076 G4HadFinalState* result = nullptr;
0077 G4int reentryCount = 0;
0078
0079 do {
0080 try {
0081
0082 result = anInteraction->ApplyYourself(thePro, *target);
0083 ++reentryCount;
0084 } catch (G4HadronicException& aR) {
0085 G4ExceptionDescription ed;
0086 aR.Report(ed);
0087 ed << "Call for " << anInteraction->GetModelName() << G4endl;
0088 ed << "Target element " << anElement->GetName() << " Z= " << target->GetZ_asInt()
0089 << " A= " << target->GetA_asInt() << G4endl;
0090 DumpState(aTrack, "ApplyYourself", ed);
0091 ed << " ApplyYourself failed" << G4endl;
0092 G4Exception("CMSSIMPInelasticProcess::PostStepDoIt", "had006", FatalException, ed);
0093 }
0094
0095
0096 result = CheckResult(thePro, *target, result);
0097 if (reentryCount > 100) {
0098 G4ExceptionDescription ed;
0099 ed << "Call for " << anInteraction->GetModelName() << G4endl;
0100 ed << "Target element " << anElement->GetName() << " Z= " << target->GetZ_asInt()
0101 << " A= " << target->GetA_asInt() << G4endl;
0102 DumpState(aTrack, "ApplyYourself", ed);
0103 ed << " ApplyYourself does not completed after 100 attempts" << G4endl;
0104 G4Exception("CMSSIMPInelasticProcess::PostStepDoIt", "had006", FatalException, ed);
0105 }
0106 } while (!result);
0107
0108
0109
0110
0111
0112 G4int nSec = result->GetNumberOfSecondaries();
0113 if (nSec > 0) {
0114 for (G4int i = 0; i < nSec; ++i) {
0115 G4DynamicParticle* dynamicParticle = result->GetSecondary(i)->GetParticle();
0116 const G4ParticleDefinition* particleDefinition = dynamicParticle->GetParticleDefinition();
0117 if (particleDefinition == G4KaonZero::Definition() || particleDefinition == G4AntiKaonZero::Definition()) {
0118 G4ParticleDefinition* newPart;
0119 if (G4UniformRand() > 0.5) {
0120 newPart = G4KaonZeroShort::Definition();
0121 } else {
0122 newPart = G4KaonZeroLong::Definition();
0123 }
0124 dynamicParticle->SetDefinition(newPart);
0125 G4ExceptionDescription ed;
0126 ed << " Hadronic model " << anInteraction->GetModelName() << G4endl;
0127 ed << " created " << particleDefinition->GetParticleName() << G4endl;
0128 ed << " -> forced to be " << newPart->GetParticleName() << G4endl;
0129 G4Exception("G4HadronicProcess::PostStepDoIt", "had007", JustWarning, ed);
0130 }
0131 }
0132 }
0133
0134 result->SetTrafoToLab(thePro.GetTrafoToLab());
0135
0136 ClearNumberOfInteractionLengthLeft();
0137
0138 FillResult(result, aTrack);
0139
0140 return theTotalResult;
0141 }