Commit 9ecebac5 authored by Andrei Gheata's avatar Andrei Gheata
Browse files

Added protection for tracks looping in field.

parent 6653bdb7
Pipeline #1060119 passed with stage
in 1 minute and 21 seconds
......@@ -134,6 +134,7 @@ private:
int fCharge = 0; /** Particle charge */
int fProcess = -1; /** Current process */
int fNsteps = 0; /** Number of steps made */
int fNintSteps = 0; /** Number of integration steps made in the current step */
int fMaxDepth = 0; /** Maximum geometry depth */
int fStage = 0; /** Simulation stage */
int fGeneration = 0; /** Track generation: 0=primary */
......@@ -306,6 +307,11 @@ public:
GEANT_FORCE_INLINE
int GetNsteps() const { return fNsteps; }
/** @brief Getter for the number of integration steps made in the current step */
VECCORE_ATT_HOST_DEVICE
GEANT_FORCE_INLINE
int GetNintSteps() const { return fNintSteps; }
VECCORE_ATT_HOST_DEVICE
GEANT_FORCE_INLINE
int GetMaxDepth() const { return fMaxDepth; }
......@@ -819,6 +825,16 @@ public:
GEANT_FORCE_INLINE
void IncrementNsteps(int nsteps = 1) { fNsteps += nsteps; }
/** @brief Setter for number of integration steps */
VECCORE_ATT_HOST_DEVICE
GEANT_FORCE_INLINE
void SetNintSteps(int nsteps) { fNintSteps = nsteps; }
/** @brief Increment number of steps */
VECCORE_ATT_HOST_DEVICE
GEANT_FORCE_INLINE
void IncrementNintSteps(int nsteps = 1) { fNintSteps += nsteps; }
/** @brief Setter for stage */
VECCORE_ATT_HOST_DEVICE
GEANT_FORCE_INLINE
......
......@@ -88,6 +88,7 @@ Track &Track::operator=(const Track &other)
fCharge = other.fCharge;
fProcess = other.fProcess;
fNsteps = other.fNsteps;
fNintSteps = other.fNintSteps;
fMaxDepth = other.fMaxDepth;
fStage = other.fStage;
fGeneration = other.fGeneration;
......@@ -149,6 +150,7 @@ void Track::Clear(const char *)
fCharge = 0;
fProcess = -1;
fNsteps = 0;
fNintSteps = 0;
fSpecies = kHadron;
fStatus = kAlive;
fMass = 0.;
......@@ -246,13 +248,13 @@ void Track::Print(const char *msg) const
{
const char *status[8] = {"alive", "killed", "inflight", "boundary", "exitSetup", "physics", "postponed", "new"};
printf("%s: evt=%d slt=%d part=%d prim=%d mth=%d gvc=%d eind=%d bind=%d chg=%d proc=%d nstp=%d spc=%d "
printf("%s: evt=%d slt=%d part=%d prim=%d mth=%d gvc=%d eind=%d bind=%d chg=%d proc=%d nstp=%d nintstp=%d spc=%d "
"status=%s mass=%g "
"xpos=%g ypos=%g zpos=%g xdir=%g ydir=%g zdir=%g mom=%g ene=%g time=%g pstp=%g stp=%g snxt=%g saf=%g nil=%g "
"ile=%g bdr=%d\n",
msg, fEvent, fEvslot, fParticle, fPrimaryIndx, fMother, fGVcode, fEindex, fBindex, fCharge, fProcess, fNsteps,
(int)fSpecies, status[int(fStatus)], fMass, fXpos, fYpos, fZpos, fXdir, fYdir, fZdir, fP, fE, fLocalTime,
fPstep, fStep, fSnext, fSafety, fNintLen, fIntLen, fBoundary);
fNintSteps, (int)fSpecies, status[int(fStatus)], fMass, fXpos, fYpos, fZpos, fXdir, fYdir, fZdir, fP, fE,
fLocalTime, fPstep, fStep, fSnext, fSafety, fNintLen, fIntLen, fBoundary);
TrackDataMgr::GetInstance()->PrintUserData(*this);
#ifndef VECCORE_CUDA
......
......@@ -505,6 +505,7 @@ void FieldPropagationHandler::PropagateInVolume(Track &track, double crtstep, Ta
// Reset relevant variables
track.SetStatus(kInFlight);
track.IncrementNintSteps();
track.IncreaseStep(crtstep);
track.DecreasePstep(crtstep);
......@@ -568,8 +569,7 @@ void FieldPropagationHandler::PropagateInVolume(TrackVec_t &tracks, const double
SOA3D<double> &DirectionOut = *(wsp->fDirectionOutp);
for (int itr = 0; itr < nTracks; ++itr) {
Track *pTrack = tracks[itr];
Track *pTrack = tracks[itr];
fltCharge[itr] = pTrack->Charge();
momentumMag[itr] = pTrack->P();
steps[itr] = stepSize[itr];
......@@ -627,6 +627,7 @@ void FieldPropagationHandler::PropagateInVolume(TrackVec_t &tracks, const double
// Update status, step and safety
track.SetStatus(kInFlight);
track.IncrementNintSteps();
track.IncreaseStep(stepSize[itr]);
track.DecreasePstep(stepSize[itr]);
......@@ -803,6 +804,7 @@ void FieldPropagationHandler::PropagateInVolume(TrackVec_t &tracks, const double
track.SetDirection(endDirVector);
track.NormalizeFast();
track.SetStatus(kInFlight);
track.IncrementNintSteps();
track.IncreaseStep(stepSize[itr]);
track.DecreasePstep(stepSize[itr]);
......
......@@ -41,10 +41,12 @@ void PreStepHandler::DoIt(Track *track, Basket &output, TaskData *td)
track->SetStatus(kInFlight);
track->SetPrePropagationDone(false);
track->SetPreStepVelocity(track->Velocity());
assert(track->PreStepVelocity() > 0);
// reset step length and energy deposit
track->SetStep(0.); // no setter for this member in Track
track->SetEdep(0.);
assert(track->PreStepVelocity() > 0);
// reset number of integration steps
track->SetNintSteps(0);
}
// Copy to output
output.AddTrack(track);
......
......@@ -64,10 +64,28 @@ int PropagationStage::CreateHandlers()
//______________________________________________________________________________
VECCORE_ATT_HOST_DEVICE
Handler *PropagationStage::Select(Track *track, TaskData *)
Handler *PropagationStage::Select(Track *track, TaskData *td)
{
// Retrieve the appropriate handler depending on the track charge
constexpr int kMaxIntSteps = 1000;
constexpr double kWarning_Energy = 100 * units::MeV;
constexpr double kImportant_Energy = 250 * units::MeV;
if (!fHasField || track->Charge() == 0) return fHandlers[0];
// Check if we reached the maximum number of integration steps for this track
if (track->GetNintSteps() > kMaxIntSteps) {
if (track->Ekin() < kImportant_Energy || track->GetNintSteps() > 10 * kMaxIntSteps) {
if (track->Ekin() > kWarning_Energy)
Warning("PropagateInVolume", "Track %d from event %d killed because it was looping in magnetic field",
track->Particle(), track->Event());
track->SetStatus(kKilled);
track->Stop();
track->SetStage(kSteppingActionsStage);
td->fNkilled++;
return nullptr;
}
}
return fHandlers[1];
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment