CustomPhysics/src/CMSSIMPInelasticProcess.cc

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   // if primary is not Alive then do nothing
0038   theTotalResult->Clear();
0039   theTotalResult->Initialize(aTrack);
0040   theTotalResult->ProposeWeight(aTrack.GetWeight());
0041   if (aTrack.GetTrackStatus() != fAlive) {
0042     return theTotalResult;
0043   }
0044
0045   // Find cross section at end of step and check if <= 0
0046   //
0047   G4DynamicParticle* aParticle = const_cast<G4DynamicParticle*>(aTrack.GetDynamicParticle());
0048
0049   // change this SIMP particle in a neutron
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   // Next check for illegal track status
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     // No warning for fStopButAlive which is a legal status here
0069     return theTotalResult;
0070   }
0071
0072   // Initialize the hadronic projectile from the track
0073   thePro.Initialise(aTrack);
0074   G4HadronicInteraction* anInteraction = GetHadronicInteractionList()[0];
0075
0076   G4HadFinalState* result = nullptr;
0077   G4int reentryCount = 0;
0078
0079   do {
0080     try {
0081       // Call the interaction
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     // Check the result for catastrophic energy non-conservation
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   // Check whether kaon0 or anti_kaon0 are present between the secondaries:
0108   // if this is the case, transform them into either kaon0S or kaon0L,
0109   // with equal, 50% probability, keeping their dynamical masses (and
0110   // the other kinematical properties).
0111   // When this happens - very rarely - a "JustWarning" exception is thrown.
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 }