GFlash/src/GflashHadronWrapperProcess.cc

0001 //
0002 // S.Y. Jun, August 2007
0003 //
0004 #include "SimG4Core/GFlash/interface/GflashHadronWrapperProcess.h"
0005
0006 #include "G4ios.hh"
0007 #include "G4GPILSelection.hh"
0008 #include "G4ProcessManager.hh"
0009 #include "G4ProcessVector.hh"
0010 #include "G4Track.hh"
0011 #include "G4VParticleChange.hh"
0012 #include "G4VProcess.hh"
0013
0014 using namespace CLHEP;
0015
0016 GflashHadronWrapperProcess::GflashHadronWrapperProcess(G4String processName)
0017     : particleChange(nullptr), pmanager(nullptr), fProcessVector(nullptr), fProcess(nullptr) {
0018   theProcessName = processName;
0019 }
0020
0021 GflashHadronWrapperProcess::~GflashHadronWrapperProcess() {}
0022
0023 G4VParticleChange *GflashHadronWrapperProcess::PostStepDoIt(const G4Track &track, const G4Step &step) {
0024   // process PostStepDoIt for the original process
0025
0026   particleChange = pRegProcess->PostStepDoIt(track, step);
0027
0028   // specific actions of the wrapper process
0029
0030   // update step/track information after PostStep of the original process is
0031   // done these steps will be repeated again without additional conflicts even
0032   // if the parameterized physics process doesn't take over the original process
0033
0034   particleChange->UpdateStepForPostStep(const_cast<G4Step *>(&step));
0035
0036   // we may not want to update G4Track during this wrapper process if
0037   // G4Track accessors are not used in the G4VFastSimulationModel model
0038   //  (const_cast<G4Step *> (&step))->UpdateTrack();
0039
0040   // update safety after each invocation of PostStepDoIts
0041   // G4double safety = std::max(endpointSafety - (endpointSafOrigin -
0042   // fPostStepPoint->GetPosition()).mag(),0.);
0043   // step.GetPostStepPoint()->SetSafety(safety);
0044
0045   // the secondaries from the original process are also created at the wrapper
0046   // process, but they must be deleted whether the parameterized physics is
0047   // an 'ExclusivelyForced' process or not.  Normally the secondaries will be
0048   // created after this wrapper process in G4SteppingManager::InvokePSDIP
0049
0050   // store the secondaries from ParticleChange to SecondaryList
0051
0052   G4TrackVector *fSecondary = (const_cast<G4Step *>(&step))->GetfSecondary();
0053   G4int nSecondarySave = fSecondary->size();
0054
0055   G4int num2ndaries = particleChange->GetNumberOfSecondaries();
0056   G4Track *tempSecondaryTrack;
0057
0058   for (G4int DSecLoop = 0; DSecLoop < num2ndaries; DSecLoop++) {
0059     tempSecondaryTrack = particleChange->GetSecondary(DSecLoop);
0060
0061     // Set the process pointer which created this track
0062     tempSecondaryTrack->SetCreatorProcess(pRegProcess);
0063
0064     // add secondaries from this wrapper process to the existing list
0065     fSecondary->push_back(tempSecondaryTrack);
0066
0067   }  // end of loop on secondary
0068
0069   // Now we can still impose conditions on the secondaries from ModelTrigger,
0070   // such as the number of secondaries produced by the original process as well
0071   // as those by other continuous processes - for the hadronic inelastic
0072   // interaction, it is always true that
0073   // particleChange->GetNumberOfSecondaries() > 0
0074
0075   // at this stage, updated step information after all processes involved
0076   // should be available
0077
0078   //  Print(step);
0079
0080   // call ModelTrigger for the second time - the primary loop of PostStepGPIL
0081   // is inside G4SteppingManager::DefinePhysicalStepLength().
0082
0083   pmanager = track.GetDefinition()->GetProcessManager();
0084
0085   G4double testGPIL = DBL_MAX;
0086   G4double fStepLength = 0.0;
0087   G4ForceCondition fForceCondition = InActivated;
0088
0089   fStepLength = step.GetStepLength();
0090
0091   fProcessVector = pmanager->GetPostStepProcessVector(typeDoIt);
0092
0093   // keep the current status of track, use fPostponeToNextEvent word for
0094   // this particular PostStep GPIL and then restore G4TrackStatus if the
0095   // paramterized physics doesn't meet trigger conditions in ModelTrigger
0096
0097   const G4TrackStatus keepStatus = track.GetTrackStatus();
0098
0099   (const_cast<G4Track *>(&track))->SetTrackStatus(fPostponeToNextEvent);
0100   int fpv_entries = fProcessVector->entries();
0101   for (G4int ipm = 0; ipm < fpv_entries; ipm++) {
0102     fProcess = (*fProcessVector)(ipm);
0103
0104     if (fProcess->GetProcessType() == fParameterisation) {
0105       // test ModelTrigger via PostStepGPIL
0106
0107       testGPIL = fProcess->PostStepGPIL(track, fStepLength, &fForceCondition);
0108
0109       // if G4FastSimulationModel:: ModelTrigger is true, then the parameterized
0110       // physics process takes over the current process
0111
0112       if (fForceCondition == ExclusivelyForced) {
0113         // clean up memory for changing the process - counter clean up for
0114         // the secondaries created by new G4Track in
0115         // G4HadronicProcess::FillTotalResult
0116         G4int nsec = particleChange->GetNumberOfSecondaries();
0117         for (G4int DSecLoop = 0; DSecLoop < nsec; DSecLoop++) {
0118           G4Track *tempSecondaryTrack = particleChange->GetSecondary(DSecLoop);
0119           delete tempSecondaryTrack;
0120         }
0121         particleChange->Clear();
0122
0123         // updating G4Step between PostStepGPIL and PostStepDoIt for the
0124         // parameterized process may not be necessary, but do it anyway
0125
0126         (const_cast<G4Step *>(&step))->SetStepLength(testGPIL);
0127         (const_cast<G4Track *>(&track))->SetStepLength(testGPIL);
0128
0129         step.GetPostStepPoint()->SetStepStatus(fExclusivelyForcedProc);
0130         ;
0131         step.GetPostStepPoint()->SetProcessDefinedStep(fProcess);
0132         step.GetPostStepPoint()->SetSafety(0.0);
0133
0134         // invoke PostStepDoIt: equivalent steps for
0135         // G4SteppingManager::InvokePSDIP
0136         particleChange = fProcess->PostStepDoIt(track, step);
0137
0138         // update PostStepPoint of Step according to ParticleChange
0139         particleChange->UpdateStepForPostStep(const_cast<G4Step *>(&step));
0140
0141         // update G4Track according to ParticleChange after each PostStepDoIt
0142         (const_cast<G4Step *>(&step))->UpdateTrack();
0143
0144         // update safety after each invocation of PostStepDoIts - acutally this
0145         // is not necessary for the parameterized physics process, but do it
0146         // anyway
0147         step.GetPostStepPoint()->SetSafety(0.0);
0148
0149         // additional nullification
0150         (const_cast<G4Track *>(&track))->SetTrackStatus(particleChange->GetTrackStatus());
0151       } else {
0152         // restore TrackStatus if fForceCondition !=  ExclusivelyForced
0153         (const_cast<G4Track *>(&track))->SetTrackStatus(keepStatus);
0154       }
0155       // assume that there is one and only one parameterized physics
0156       break;
0157     }
0158   }
0159
0160   // remove secondaries of this wrapper process that were added to the secondary
0161   // list since they will be added in the normal stepping procedure after
0162   // this->PostStepDoIt in G4SteppingManager::InvokePSDIP
0163
0164   // move the iterator to the (nSecondarySave+1)th element in the secondary list
0165   G4TrackVector::iterator itv = fSecondary->begin();
0166   itv += nSecondarySave;
0167
0168   // delete next num2ndaries tracks from the secondary list
0169   fSecondary->erase(itv, itv + num2ndaries);
0170
0171   // end of specific actions of this wrapper process
0172
0173   return particleChange;
0174 }
0175
0176 void GflashHadronWrapperProcess::Print(const G4Step &step) {
0177   G4cout << " GflashHadronWrapperProcess ProcessName, PreStepPosition, preStepPoint KE, PostStepPoint KE, DeltaEnergy "
0178             "Nsec \n "
0179          << step.GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName() << " "
0180          << step.GetPostStepPoint()->GetPosition() << " " << step.GetPreStepPoint()->GetKineticEnergy() / CLHEP::GeV
0181          << " " << step.GetPostStepPoint()->GetKineticEnergy() / CLHEP::GeV << " " << step.GetDeltaEnergy() / CLHEP::GeV
0182          << " " << particleChange->GetNumberOfSecondaries() << G4endl;
0183 }