Find nucleus involved in nCapture process in Geant4

Question

I'm trying to determine the target nucleus in an (n,gamma) reaction in Geant4. I have been unable to extract this information. The only place I have found it to be stored is in G4IsoParticleChange which is created by the G4HadronicProcess if Isotope Counting is enabled.

Unfortunately this information is lost (not in a nice way, either, but memory leak style) every time the process is invoked. I cannot find a user hook-in to ask for this information in between particle creations. G4UserStackingAction is not sufficient as far as I can tell because the secondary particles are all created for a step before they are stacked (* though writing that last sentence has given me an idea).

Could anyone help me determine the nucleus which captures the neutron in the (n,gamma) reaction? Is there an easier way to get this information?

Thanks

P.S. Since a neutron can only be absorbed once within a step, would it be safe to just wait until the secondaries are stacked to obtain the IsoParticleChange info or am I risking a memory leak?

Edit to be more clear:

I'm asking if someone knows how to retrieve the nucleus which was the target in the nCapture process in Geant4. It is clear from the source that a memory leak will occur whenever the G4HadronicProcess is invoked if the G4IsoParticleChange information is not retrieved. There appears to be no user hook in point which is appropriate to grab this without missing some information and causing said memory leak, and yet it is possible to turn this information storage on. I am wondering what the correct way to grab this information is, or if there is a better way to obtain this type of info about the target nucleus.

I have previously gotten some help on SO from people who are familiar with Geant. I have not been able to get a response from the Slac Geant4 forum to give me access to post there. The forum doesn't seem terribly active, anyway.

Thanks

Follow-up:

In case someone stumbles upon this, the answer I got from one of the authors of the G4HadronicProcess class was "This part of the implementation has not been maintained for many years" and "Getting the nucleus involved in the interaction is very difficult if not impossible without editing the source code." So I'm in the process of setting up a new Geant workspace where I can do that. FYI.

Answer 1

I finally found an effective solution to that problem: I search for a nucleus in the secondary track vector of a step ending by a neutron capture. This nucleus is the one which captured the neutron and which recoil is tracked by Geant4. Do not forget to subtract 1 neutron to this nucleus to get what you want (for a capture on hydrogen, you will get deuterium with this method).

So in my SteppingAction, I add:

// neutron capture
if ( aStep->GetPostStepPoint()->GetProcessDefinedStep()
  && aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessType() == fHadronic // see G4ProcessType.hh
  && aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessSubType() == fCapture // see G4HadronicProcessType.hh
) {
  if ( aStep->GetSecondary() != 0
    && aStep->GetSecondary()->size() != 0
  ) {
    std::vector<G4Track*>::const_iterator it;
    for (it=aStep->GetSecondary()->begin(); it!=aStep->GetSecondary()->end(); it++) {
    if ( !(*it)->GetCreatorProcess()
      ||  (*it)->GetCreatorProcess()->GetProcessSubType() != fCapture // see G4HadronicProcessType.hh
      ||  (*it)->GetCreatorProcess()->GetProcessType() != fHadronic // see G4ProcessType.hh
      || !(*it)->GetDynamicParticle()
      ||  (*it)->GetDynamicParticle()->GetParticleDefinition()->GetAtomicNumber() == 0 // keep only nucleus
      ||  (*it)->GetDynamicParticle()->GetParticleDefinition() == G4Neutron::NeutronDefinition() // but not the neutron (perhaps antiparticles could also be checked)
    ) { continue; }
      myEventAction->TreatNeutronCapture(*it);
      break;
    }
    if (it == aStep->GetSecondary()->end()) {
      G4cerr << "\n#### END OF SECONDARY VECTOR REACHED AFTER NEUTRON CAPTURE ! ###\n";
      myEventAction->TreatNeutronCapture(0);
    }
  } else { myEventAction->TreatNeutronCapture(0); }
}

As you can see, I defer the actual stepping treatment to methods of EventAction. This is my TreatNeutronCapture method :

void MyExperimentEventAction::TreatNeutronCapture(const G4Track* track)
{
  myParticle = myMC->GetParticle(ParticleMap[track->GetParentID()]);
  if (track == 0) { myParticle->SetFinalProcess(-1); } // if no nucleus is found in secondaries
  else { myParticle->SetFinalProcess(track->GetDynamicParticle()->GetPDGcode()); }
}

myParticle and myMC refers to my own data classes.

Answer 2

Well, I haven't looked in detail to see if it actually meets you needs, but it seems that many problems related to keeping track of parents and products and the like can be solved in a manner analogous to that described in this Tip from the forum.

Aha! The tip on examining secondaries may be better because it shows how to select on the ending state of the parent particle and extract some information from the process that ended the step.

You might start there.