Problem 726

Summary: Possible bug in G4LCapture
Product: Geant4 Reporter: seligman
Component: processes/hadronic/modelsAssignee: dennis.herbert.wright
Status: RESOLVED FIXED    
Severity: normal CC: dennis.herbert.wright
Priority: P2    
Version: 7.0   
Hardware: PC   
OS: Linux   

Description seligman 2005-02-24 14:18:46 CET 
While doing detailed calibration studies of energy deposited in ATLAS, we see
signs of a possible bug in G4LCapture.

We ran a study with 200-GeV pions in the ATLAS detector.  It occurred in about
1% of those events, and only in 1-2 G4Steps in such events.  Hence, this issue
is quite rare.

What appears to happen is that in a low-energy neutron capture, only 1-2 gammas
with a few MeV is are supposed to be emitted.  But in the rare G4Steps, the
energy of those gammas is on the order of several GeV.

When I look at the code in G4LCapture, it's not clear how this is supposed to
happen.  The code seems designed to emit gammas with energies on the order of a
few MeV.  Something anomalous is going on.

I will append more detailed diagnostic information in a follow-up to this
initial bug report.
Comment 1 seligman 2005-02-24 14:21:59 CET 
This is a diagnostic that we generated from information available from G4Step,
the G4SteppingManager, and G4VProcess.  The sense of the error messages is that
if we detect "neative invisible energy" (a neutrino with a few MeV comes in,
several GeV come out), we print out detailed information about all the
secondaries in the G4Step.

SimulationEnergies: negative invisible energy (3) --
  event=1  phys volume=LAr::EMEC::FrontOuterLongBar
  particle name=neutron mass=939.5656 material=Gten status=2
midStepGlobalTime=22.4275
  EkinPreStep=4.8229 EkinPostStep=0.0000 dEkinStep=4.8229
  dEStepVisible=0.0000 dEkinNonDeposited=4.8229 EtotalPostStep=939.5656
  secondaryParticleName=gamma PDGcode=22 totalEofSecondary=690.5261
kinEofSecondary=690.5261 result.energy[kInvisible0]=-4389.3471
 processDefinedStepName=LCapture ProcessType=NotDefined
  fN2ndariesAtRestDoIt=0 fN2ndariesAlongStepDoIt=0 fN2ndariesPostStepDoIt=2
tN2ndariesTot=2
  loopStart=0 loopEnd=2 lp1=1 globalTimeOfSecondary=22.4642
================= Detailed report ================
 loopStart=0 loopEnd=2
 i=0: name=gamma PDG=22
   totalE=3703.6439 kinE=3703.6439 time=22.4642 vertexKinE=0.0000
   creatorProcess=LCapture
 i=1: name=gamma PDG=22
   totalE=690.5261 kinE=690.5261 time=22.4642 vertexKinE=0.0000
   creatorProcess=LCapture
================= End of detailed report ================

 G4Transportation is killing track that is looping or stuck
   This track has 0.33447043 MeV energy.

SimulationEnergies: negative invisible energy (3) --
  event=1  phys volume=LAr::EMB::ThinAbs::CornerDownFold
  particle name=neutron mass=939.5656 material=Thinabs status=2
midStepGlobalTime=9.0489
  EkinPreStep=10.3791 EkinPostStep=0.0000 dEkinStep=10.3791
  dEStepVisible=0.0000 dEkinNonDeposited=10.3791 EtotalPostStep=939.5656
  secondaryParticleName=gamma PDGcode=22 totalEofSecondary=8388.1634
kinEofSecondary=8388.1634 result.energy[kInvisible0]=-8377.7844
 processDefinedStepName=LCapture ProcessType=NotDefined
  fN2ndariesAtRestDoIt=0 fN2ndariesAlongStepDoIt=0 fN2ndariesPostStepDoIt=2
tN2ndariesTot=2
  loopStart=0 loopEnd=2 lp1=0 globalTimeOfSecondary=9.1132

SimulationEnergies: negative invisible energy (3) --
  event=1  phys volume=LAr::EMB::ThinAbs::CornerDownFold
  particle name=neutron mass=939.5656 material=Thinabs status=2
midStepGlobalTime=9.0489
  EkinPreStep=10.3791 EkinPostStep=0.0000 dEkinStep=10.3791
  dEStepVisible=0.0000 dEkinNonDeposited=10.3791 EtotalPostStep=939.5656
  secondaryParticleName=gamma PDGcode=22 totalEofSecondary=4226.4416
kinEofSecondary=4226.4416 result.energy[kInvisible0]=-12604.2260
 processDefinedStepName=LCapture ProcessType=NotDefined
  fN2ndariesAtRestDoIt=0 fN2ndariesAlongStepDoIt=0 fN2ndariesPostStepDoIt=2
tN2ndariesTot=2
  loopStart=0 loopEnd=2 lp1=1 globalTimeOfSecondary=9.1132
================= Detailed report ================
 loopStart=0 loopEnd=2
 i=0: name=gamma PDG=22
   totalE=8388.1634 kinE=8388.1634 time=9.1132 vertexKinE=0.0000
   creatorProcess=LCapture
 i=1: name=gamma PDG=22
   totalE=4226.4416 kinE=4226.4416 time=9.1132 vertexKinE=0.0000
   creatorProcess=LCapture
================= End of detailed report ================



SimulationEnergies: negative invisible energy (3) --
  event=1  phys volume=LAr::EMEC::BackOuterLongBar
  particle name=neutron mass=939.5656 material=Gten status=2
midStepGlobalTime=16.7194
  EkinPreStep=30.4551 EkinPostStep=0.0000 dEkinStep=30.4551
  dEStepVisible=0.0000 dEkinNonDeposited=30.4551 EtotalPostStep=939.5656
  secondaryParticleName=gamma PDGcode=22 totalEofSecondary=14398.2297
kinEofSecondary=14398.2297 result.energy[kInvisible0]=-14367.7746
 processDefinedStepName=LCapture ProcessType=NotDefined
  fN2ndariesAtRestDoIt=0 fN2ndariesAlongStepDoIt=0 fN2ndariesPostStepDoIt=2
tN2ndariesTot=2
  loopStart=2 loopEnd=4 lp1=2 globalTimeOfSecondary=16.7809

SimulationEnergies: negative invisible energy (3) --
  event=1  phys volume=LAr::EMEC::BackOuterLongBar
  particle name=neutron mass=939.5656 material=Gten status=2
midStepGlobalTime=16.7194
  EkinPreStep=30.4551 EkinPostStep=0.0000 dEkinStep=30.4551
  dEStepVisible=0.0000 dEkinNonDeposited=30.4551 EtotalPostStep=939.5656
  secondaryParticleName=gamma PDGcode=22 totalEofSecondary=124.1505
kinEofSecondary=124.1505 result.energy[kInvisible0]=-14491.9252
 processDefinedStepName=LCapture ProcessType=NotDefined
  fN2ndariesAtRestDoIt=0 fN2ndariesAlongStepDoIt=0 fN2ndariesPostStepDoIt=2
tN2ndariesTot=2
  loopStart=2 loopEnd=4 lp1=3 globalTimeOfSecondary=16.7809
================= Detailed report ================
 loopStart=2 loopEnd=4
 i=2: name=gamma PDG=22
   totalE=14398.2297 kinE=14398.2297 time=16.7809 vertexKinE=0.0000
   creatorProcess=LCapture
 i=3: name=gamma PDG=22
   totalE=124.1505 kinE=124.1505 time=16.7809 vertexKinE=0.0000
   creatorProcess=LCapture
================= End of detailed report ================
Comment 2 seligman 2005-02-24 14:25:59 CET 
For technical reasons associated with the ATLAS Athena framework, I could not
place debug statements within G4LCapture.  This is an example of a revised
G4LCapture which I believe would generate diagnostics if the problem were to
repeat itself in a test simulation.

//
// ********************************************************************
// * DISCLAIMER                                                       *
// *                                                                  *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in:                      *
// *   http://cern.ch/geant4/license                                  *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.                                                             *
// *                                                                  *
// * This  code  implementation is the  intellectual property  of the *
// * GEANT4 collaboration.                                            *
// * By copying,  distributing  or modifying the Program (or any work *
// * based  on  the Program)  you indicate  your  acceptance of  this *
// * statement, and all its terms.                                    *
// ********************************************************************
//
//
// $Id: G4LCapture.cc,v 1.10 2004/12/07 13:49:06 gunter Exp $
// GEANT4 tag $Name: geant4-07-00-cand-03 $
//
//
// G4 Model: Low-energy Neutron Capture
// F.W. Jones, TRIUMF, 03-DEC-96
//
// This is a prototype of a low-energy neutron capture process.
// Currently it is based on the GHEISHA routine CAPTUR,
// and conforms fairly closely to the original Fortran.
//
// HPW Capture using models now. the code comes from the
// original G4LCapture class.
//
// 25-JUN-98 FWJ: replaced missing Initialize for ParticleChange.
//

#include "globals.hh"
#include "G4LCapture.hh"
#include "Randomize.hh"

G4LCapture::G4LCapture() :
   G4HadronicInteraction()
{
   SetMinEnergy( 0.0*GeV );
   SetMaxEnergy( DBL_MAX );
}

G4LCapture::~G4LCapture()
{
   theParticleChange.Clear();
}

G4HadFinalState*
G4LCapture::ApplyYourself(const G4HadProjectile & aTrack, G4Nucleus& targetNucleus)
{

   theParticleChange.Clear();
   theParticleChange.SetStatusChange(stopAndKill);

   G4double N = targetNucleus.GetN();
   G4double Z = targetNucleus.GetZ();

   const G4LorentzVector theMom = aTrack.Get4Momentum();
   G4double P = theMom.vect().mag()/GeV;
   G4double Px = theMom.vect().x();
   G4double Py = theMom.vect().y();
   G4double Pz = theMom.vect().z();
   G4double E = theMom.e()/GeV;
   G4double mass = aTrack.GetDefinition()->GetPDGMass()/GeV;
   G4double Q = aTrack.GetDefinition()->GetPDGCharge();
   if (verboseLevel > 1) {
      G4cout << "G4LCapture:ApplyYourself: incident particle:" << G4endl;
      G4cout << "P      " << P << " GeV/c" << G4endl;
      G4cout << "Px     " << Px << " GeV/c" << G4endl;
      G4cout << "Py     " << Py << " GeV/c" << G4endl;
      G4cout << "Pz     " << Pz << " GeV/c" << G4endl;
      G4cout << "E      " << E << " GeV" << G4endl;
      G4cout << "mass   " << mass << " GeV" << G4endl;
      G4cout << "charge " << Q << G4endl;
   }
// GHEISHA ADD operation to get total energy, mass, charge:
   if (verboseLevel > 1) {
      G4cout << "G4LCapture:ApplyYourself: material:" << G4endl;
      G4cout << "A      " << N << G4endl;
      G4cout << "Z   " << Z  << G4endl;
      G4cout << "atomic mass " <<
        Atomas(N, Z) << "GeV" << G4endl;
   }
   E = E + Atomas(N, Z);
   G4double E02 = E*E - P*P;
   G4double E0 = std::sqrt(std::abs(E02));
   if (E02 < 0) E0 = -E0;
   Q = Q + Z;
   if (verboseLevel > 1) {
      G4cout << "G4LCapture:ApplyYourself: total:" << G4endl;
      G4cout << "E      " << E << " GeV" << G4endl;
      G4cout << "mass   " << mass << " GeV" << G4endl;
      G4cout << "charge " << Q << G4endl;
   }
   Px = -Px;
   Py = -Py;
   Pz = -Pz;

// Make a gamma...

   G4double ran = G4RandGauss::shoot();
   G4double p = 0.0065 + ran*0.0010;

   G4double ran1 = G4UniformRand();
   G4double ran2 = G4UniformRand();
   G4double cost = -1. + 2.*ran1;
   G4double sint = std::sqrt(std::abs(1. - cost*cost));
   G4double phi = ran2*twopi;

   G4double px = p*sint*std::sin(phi);
   G4double py = p*sint*std::cos(phi);
   G4double pz = p*cost;

   double savepx = px;
   double savepy = py;
   double savepz = pz;

   G4double e = p;
   G4double e0 = 0.;

   G4double a = px*Px + py*Py + pz*Pz;
   a = (a/(E + E0) - e)/E0;

   px = px + a*Px;
   py = py + a*Py;
   pz = pz + a*Pz;

   G4DynamicParticle* aGamma;
   aGamma = new G4DynamicParticle(G4Gamma::GammaDefinition(),
                                  G4ThreeVector(px*GeV, py*GeV, pz*GeV));
   theParticleChange.AddSecondary(aGamma);

   bool problem = false;
   double egamma = px*px + py*py + pz*pz;
   if ( egamma > 3.*GeV || p < 0. )
     {
       problem = true;
       std::cout << std::endl << "============ G4LCapture Debug
=================" << std::endl;
      G4cout << "G4LCapture:ApplyYourself: incident particle:" << G4endl;
      G4cout << "P      " << P << " GeV/c" << G4endl;
      G4cout << "Px     " << Px << " GeV/c" << G4endl;
      G4cout << "Py     " << Py << " GeV/c" << G4endl;
      G4cout << "Pz     " << Pz << " GeV/c" << G4endl;
      G4cout << "E      " << E << " GeV" << G4endl;
      G4cout << "mass   " << mass << " GeV" << G4endl;
      G4cout << "charge " << Q << G4endl;
      G4cout << "G4LCapture:ApplyYourself: material:" << G4endl;
      G4cout << "A      " << N << G4endl;
      G4cout << "Z   " << Z  << G4endl;
      G4cout << "atomic mass " <<
        Atomas(N, Z) << "GeV" << G4endl;
      G4cout << "G4LCapture:ApplyYourself: total:" << G4endl;
      G4cout << "E      " << E << " GeV" << G4endl;
      G4cout << "E0   " << E0 << " GeV" << G4endl;
      G4cout << "charge " << Q << G4endl;
      std::cout << "additional debug:" << std::endl;
      std::cout << " p=" << p << std::endl
                << " save (px,py,pz)=(" << savepx << "," << savepy << "," <<
savepz << ")" << std::endl;
      std::cout << " a=" << a << std::endl
                << " (px,py,pz)=(" << px << "," << py << "," << pz << ")" <<
std::endl;
       std::cout << "-----------------------------------------------" << std::endl;
     }


// Make another gamma if there is sufficient energy left over...

   G4double xp = 0.008 - p;
   if (xp <= 0.) return &theParticleChange;

   ran1 = G4UniformRand();
   ran2 = G4UniformRand();
   cost = -1. + 2.*ran1;
   sint = std::sqrt(std::abs(1. - cost*cost));
   phi = ran2*twopi;

   px = xp*sint*std::sin(phi);
   py = xp*sint*std::cos(phi);
   pz = xp*cost;

   savepx = px;
   savepy = py;
   savepz = pz;

   e = xp;
   e0 = 0.;

   a = px*Px + py*Py + pz*Pz;
   a = (a/(E + E0) - e)/E0;

   px = px + a*Px;
   py = py + a*Py;
   pz = pz + a*Pz;

   aGamma = new G4DynamicParticle(G4Gamma::GammaDefinition(),
                                  G4ThreeVector(px*GeV, py*GeV, pz*GeV));
   theParticleChange.AddSecondary(aGamma);

   if ( problem )
     {
      std::cout << "additional debug, second gamma:" << std::endl;
      std::cout << " xp=" << xp << std::endl
                << " save (px,py,pz)=(" << savepx << "," << savepy << "," <<
savepz << ")" << std::e
ndl;
      std::cout << " a=" << a << std::endl
                << " (px,py,pz)=(" << px << "," << py << "," << pz << ")" <<
std::endl;
      std::cout << "===============================================" <<
std::endl << std::endl;
     }

   return &theParticleChange;
}
Comment 3 dennis.herbert.wright 2005-03-15 18:18:59 CET 
The problem is due to incorrect units in G4LCapture::ApplyYourself.  The 7th,
8th and 9th lines of that method,
   G4double Px = theMom.vect().x();
   G4double Py = theMom.vect().y();
   G4double Pz = theMom.vect().z();
should be replaced with
   G4double Px = theMom.vect().x()/GeV;
   G4double Py = theMom.vect().y()/GeV;
   G4double Pz = theMom.vect().z()/GeV;

Also, although it is not required for the fix, it is easy to put in the correct
physics for neutron capture on the proton:
beginning at the line:

 // Make a gamma...

replace

   G4double ran = G4RandGauss::shoot();
   G4double p = 0.0065 + ran*0.0010;

with

   G4double p;
   if (Z == 1) {  // special case for hydrogen
     p = 0.0022;
   } else {
     G4double ran = G4RandGauss::shoot();
     p = 0.0065 + ran*0.0010;
   }

and put
   if ( Z > 1 ) {

   }

around the code which makes a second gamma.
Comment 4 Hans-Peter.Wellisch 2005-07-05 07:01:59 CEST 
thank you for reporting this.