Decay table problem with G4

[sawenzel]

I am running into fatal core dumps when running with G4:

-------- EEEE ------- G4Exception-START -------- EEEE ------- *** G4Exception : DECAY003 issued by : G4Decay::DoIt Can not determine decay channel for sigma+ mass of dynamic particle: 1.05831 (GEV) dacay table has 2 entries 0: BR 0.516, IsOK? 0, --> proton + pi0 1: BR 0.483, IsOK? 0, --> neutron + pi+

Fatal Exception core dump ***

Are we missing any trivial configuration? How can this be fixed?

[sawenzel]

@ihrivnac, @preghenella , @shahor02 : Any idea?

[ihrivnac]

What's your Geant4 version?

[sawenzel]

It is v10.3.3 (as specified in defaults-o2.sh).

[ihrivnac]

It looks that the problem is in the " mass of dynamic particle: 1.05831 (GEV)" which is smaller than the PDG mass, which is for sigma+ 1.18937*GeV. Could you investigate if it's a primary or a secondary? In the first case, the problem is in the event generator, and in the second case in Geant4 physics process.

[sawenzel]

Here is the track as seen in gdb:

p *track
$2 = {fCurrentStepNumber = 254, fPosition = {dx = -0.61883177255248523, dy = -8.4389969345702642, dz = -28145.101366639814, static tolerance = 2.22045e-14}, fGlobalTime = 93.881970651531319, 
  fLocalTime = 93.881970651531319, fTrackLength = 19098.602738877464, fParentID = 0, fTrackID = 1689, fVelocity = 299.79245800000001, fpTouchable = {fObj = 0xa59a470}, fpNextTouchable = {fObj = 0xa59a470}, 
  fpOriginTouchable = {fObj = 0x10db6bd0}, fpDynamicParticle = 0x10a33c80, fTrackStatus = fAlive, fBelowThreshold = false, fGoodForTracking = false, fStepLength = 9046.5040046367649, fWeight = 1, fpStep = 0xa531160, 
  fVtxPosition = {dx = 0, dy = 0, dz = 0, static tolerance = 2.22045e-14}, fVtxMomentumDirection = {dx = -1.8036656520484252e-05, dy = -0.00022300838818432419, dz = -0.9999999749709686, static tolerance = 2.22045e-14}, 
  fVtxKineticEnergy = 1733876.0118748124, fpLVAtVertex = 0xbac1170, fpCreatorProcess = 0x0, fCreatorModelIndex = -1, fpUserInformation = 0x10ae27e0, prev_mat = 0x0, groupvel = 0x0, prev_velocity = 0, prev_momentum = 0, 
  is_OpticalPhoton = false, static velTable = 0x1087f010, useGivenVelocity = false, fpAuxiliaryTrackInformationMap = 0x0}

Maybe you can spot if primary or not; Otherwise I need to debug a bit deeper.

[ihrivnac]

The "fParentID = 0" means, that it is primary.

[sawenzel]

This is the (unmodified) primary from the event (as stored inside stack):

p *(mPrimaryParticles._M_impl._M_start + 1689)
$14 = {<TObject> = {_vptr.TObject = 0x7fffed37f008 <vtable for TParticle+16>, fUniqueID = 0, fBits = 33554432, static fgDtorOnly = 0, static fgObjectStat = false, static fgIsA = {_M_b = {
        _M_p = 0x145f990}}}, <TAttLine> = {_vptr.TAttLine = 0x7fffed37f1f0 <vtable for TParticle+504>, fLineColor = 1, fLineStyle = 1, fLineWidth = 1, static fgIsA = {_M_b = {_M_p = 0x1f46b10}}}, <TAtt3D> = {
    _vptr.TAtt3D = 0x7fffed37f280 <vtable for TParticle+648>, static fgIsA = {_M_b = {_M_p = 0x175ff50}}}, fPdgCode = -211, fStatusCode = 1689, fMother = {-1, 0}, fDaughter = {-1, -1}, fWeight = 1, fCalcMass = 0.13957, 
  fPx = 0.19439059495925903, fPy = 0.14248257875442505, fPz = -290.89974975585938, fE = 290.89987182617188, fVx = 0, fVy = 0, fVz = 0, fVt = 0, fPolarTheta = -99, fPolarPhi = -99, fParticlePDG = 0xa31d050, 
  static fgIsA = {_M_b = {_M_p = 0x3f4c8c0}}}

[ihrivnac]

This does not look to be the same primary, as this one has fPdgCode = -211, what means that it is pi-, while the particle which caused exception is sigma+.

[sawenzel]

Right, G4 seems to shift the track indices by 1 which is why I picked the wrong one. This should be the correct one:

 p (*(o2::Data::Stack*)(this->fMC->fStack)).mCurrentParticle
$3 = {<TObject> = {_vptr.TObject = 0x7fffed37f008 <vtable for TParticle+16>, fUniqueID = 0, fBits = 50331648, static fgDtorOnly = 0, static fgObjectStat = false, static fgIsA = {_M_b = {_M_p = 0x145f940}}}, <TAttLine> = {
    _vptr.TAttLine = 0x7fffed37f1f0 <vtable for TParticle+504>, fLineColor = 1, fLineStyle = 1, fLineWidth = 1, static fgIsA = {_M_b = {_M_p = 0x1f46960}}}, <TAtt3D> = {
    _vptr.TAtt3D = 0x7fffed37f280 <vtable for TParticle+648>, static fgIsA = {_M_b = {_M_p = 0x175fec0}}}, fPdgCode = 3222, fStatusCode = 1688, fMother = {-1, 0}, fDaughter = {-1, -1}, fWeight = 1, fCalcMass = 1.18937, 
  fPx = -0.03129240870475769, fPy = -0.38690483570098877, fPz = -1734.9339599609375, fE = 1734.934326171875, fVx = 0, fVy = 0, fVz = 0, fVt = 0, fPolarTheta = -99, fPolarPhi = -99, fParticlePDG = 0xa333830, 
  static fgIsA = {_M_b = {_M_p = 0x3f4c8d0}}}

The mass of this particle seems to be well 1.18 GeV.

[sawenzel]

I am also adding the instances of G4 dynamic particle (at the moment of crash) and the original primary:

p *track->fpDynamicParticle
$6 = {theMomentumDirection = {dx = -9.9320526059858683e-05, dy = 0.0016063573805523157, dz = -0.99999870487486076, static tolerance = 2.22045e-14}, theParticleDefinition = 0xcf14110, thePolarization = {dx = 0, dy = 0, 
    dz = 0, static tolerance = 2.22045e-14}, theKineticEnergy = 1733873.8702146623, theProperTime = 0.057268237925851805, theDynamicalMass = 1058.3142970626388, theDynamicalCharge = 1, theDynamicalSpin = 0.5, 
  theDynamicalMagneticMoment = 7.7487251332839014e-11, theElectronOccupancy = 0x0, thePreAssignedDecayProducts = 0x0, thePreAssignedDecayTime = -1, verboseLevel = 1, primaryParticle = 0x10896800, thePDGcode = 0}

(gdb) p *track->fpDynamicParticle->primaryParticle
$7 = {_vptr.G4PrimaryParticle = 0x7fffd9b29968 <vtable for G4PrimaryParticle+16>, PDGcode = 3222, G4code = 0xcf14110, direction = {dx = -1.8036656520484252e-05, dy = -0.00022300838818432419, dz = -0.9999999749709686, 
    static tolerance = 2.22045e-14}, kinE = 1733876.0118748124, nextParticle = 0x10896890, daughterParticle = 0x0, trackID = 1689, mass = 1058.3142970626388, charge = 1, polX = 0, polY = 0, polZ = 0, Weight0 = 1, 
  properTime = 0, userInfo = 0x0}

So, could it be that the mass gets incorrectly converted during setup already?

[ihrivnac]

We do not pass static mass (and other static properties) from TParticle to G4 particle; we just get the pdg code and let G4 find G4ParticleDefinition for this PDG code. In G4, the mass defined for sigma+ (3222) is well 1.18937, so a question is where this primary gets 1058. Could you extract the properties of the "G4code = 0xcf14110" object, which should be the G4ParticleDefinition for sigma+? Also check, if there are not warnings about not found particle definitions: "G4ParticleTable::FindParticle() for particle with pdgEncoding= ... failed"

[sawenzel]

Here is the G4code object:

p *track->fpDynamicParticle->primaryParticle->G4code
$8 = {_vptr.G4ParticleDefinition = 0x7fffd9b296f0 <vtable for G4ParticleDefinition+16>, theProcessManagerShadow = 0xd3e4170, g4particleDefinitionInstanceID = 428, static subInstanceManager = {
    static slavetotalspace = 513, static offset = 0xd02f3c0, totalobj = 464, mutex = 0}, theQuarkContent = {0, 2, 1, 0, 0, 0}, theAntiQuarkContent = {0, 0, 0, 0, 0, 0}, 
  theParticleName = {<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >> = {static npos = 18446744073709551615, 
      _M_dataplus = {<std::allocator<char>> = {<__gnu_cxx::new_allocator<char>> = {<No data fields>}, <No data fields>}, _M_p = 0xcf14168 "sigma+"}, _M_string_length = 6, {
        _M_local_buf = "sigma+\000\000\000\000\000\000\000\000\000", _M_allocated_capacity = 47697447315827}}, <No data fields>}, thePDGMass = 1189.3700000000001, thePDGWidth = 8.2089999999999995e-12, thePDGCharge = 1, 
  thePDGiSpin = 1, thePDGSpin = 0.5, thePDGiParity = 1, thePDGiConjugation = 0, thePDGiGParity = 0, thePDGiIsospin = 2, thePDGiIsospin3 = 2, thePDGIsospin = 1, thePDGIsospin3 = 1, 
  thePDGMagneticMoment = 7.7487251332839014e-11, theLeptonNumber = 0, theBaryonNumber = 1, theParticleType = {<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >> = {
      static npos = 18446744073709551615, _M_dataplus = {<std::allocator<char>> = {<__gnu_cxx::new_allocator<char>> = {<No data fields>}, <No data fields>}, _M_p = 0xcf141e8 "baryon"}, _M_string_length = 6, {
        _M_local_buf = "baryon\000\000\000\000\000\000\000\000\000", _M_allocated_capacity = 121425057964386}}, <No data fields>}, 
  theParticleSubType = {<std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >> = {static npos = 18446744073709551615, 
      _M_dataplus = {<std::allocator<char>> = {<__gnu_cxx::new_allocator<char>> = {<No data fields>}, <No data fields>}, _M_p = 0xcf14208 "sigma"}, _M_string_length = 5, {
        _M_local_buf = "sigma\000\000\000\000\000\000\000\000\000\000", _M_allocated_capacity = 418447321459}}, <No data fields>}, thePDGEncoding = 3222, theAntiPDGEncoding = -3222, fShortLivedFlag = false, 
  thePDGStable = false, thePDGLifeTime = 0.080180000000000001, theDecayTable = 0xcf142f0, theParticleTable = 0x7fffd9b2f6c0 <G4ParticleTable::GetParticleTable()::theParticleTable>, theAtomicNumber = 0, theAtomicMass = 0, 
  verboseLevel = 1, fApplyCutsFlag = false, isGeneralIon = false}

I did not find "pdgEncoding" warnings.

[ihrivnac]

This looks correct. Do you get this fatal with O2 or AliRoot?

[sawenzel]

Its with O2. Do you have any hint where to set a breakpoint in order to find out if something goes wrong with the G4DynamicParticles (or similar) creation?

[ihrivnac]

I am building a "reference" version (which means with alidist defaults, but Geant4 MT on), when done I will try to reproduce the problem. To investigate, where the mass gets corrupted, I suggest to add printing of the mass at stepping level when a track has pdg = 3222. How long does the program runs to get the crash?

[ihrivnac]

My installation has finished, I don't get this fatal when running with: o2sim --modules EMCAL TPC TOF TRD --mcEngine TGeant4 --nEvents 100 How do you run your code?

[sawenzel]

I will set up a reproducible run for you. I see my crash with a Hijing event which I might need to pass to you.

[sawenzel]

I am seeing the error (both on linux + mac) when running

o2sim -m TRD ITS TPC PIPE -e TGeant4 -n 1 -g extkin --extKinFile Kinematics_Hijing_small.root

I will send you the kinematics file by email (cannot attach here unfortunately).

[ihrivnac]

Thank you. I can reproduce the problem, but it's not easy to debug, as when using the standard O2 build procedure, each modification of geant4_vmc (defined as a development package), triggers re-build of FairRoot and O2 ...

[sawenzel]

You can avoid rebuilding everything by checking out geant4_vmc as a development package and rebuild it using make install inside sw/BUILD/geant4_vmc/.../. This will just update the installed libraries.

[ihrivnac]

This is what I did first, but modifications seemed not to be taken into account. I will try it out with the next iteration.

[ihrivnac]

In difference from what I thought (and wrote above), the mass of the primary particle is computed from the input momentum and energy [ a square root from (EE - pp)], as we use the ctor which takes all these values: G4PrimaryParticle(const G4ParticleDefinition* Gcode, G4double px,G4double py,G4double pz,G4double E);

In this wrong case, the input values (in MeV) of the primary passed from the generator are: momentum x: -31.2924
momentum y: -386.905 momentum z: -1.73493e+06 energy: 1.73493e+06 so the result looks to be affected by the precision of the computation.

This code in both Geant4 VMC and Geant4 seems to be unchanged since the versions which we use in productions, where such fatal was not reported, Did you make some changes in O2 in event generator (eg. passing from double to float) for optimization reasons ?

[sawenzel]

At this level everything is in double precision. I analysed the kinematics file, and it turns out that the sigma particle is already stored like this and has the discrepancy of calculated vs tabulated mass. The file is obtained straight from the Hijing generator and AliRoot kinematics output. Either this output process does something to the precision or we were just lucky not to see it in AliRoot.

What options do we have? a) filter out particles that have a mass descrepancy? b) rescale momentum, so that it matches? ...

[pzhristov]

I think the right solution is to fix the generator. Lets have a look why we get particles off the mass shell.

[ihrivnac]

Yes, I think the best would be to understand the problem with the generator first. If needed, we can put a temporary work-around, put a check of the mass in Geant4 VMC, print a warning here, so that such cases are still detected, and then pass only momentum to G4 primary constructor.

[sawenzel]

Dear All. Thanks for your efforts and input here. I will study why we get these particles from the generator. For the moment, we don't need to do something urgently (I can simply take out the problematic particles).

Maybe we can add some special treatment in G4 VMC later.

[sawenzel]

Another question: Could it be that this effect is not seen in AliRoot because there we attach external decayers?

[ihrivnac]

The external decayer is applied only for particles which do not have defined a decay table in Geant4, so this exception should happen in AliRoot, too. Are there modifications in the event generator wrt AliRoot generator?