Get proper hit/MCParticle relations from Geant4Output2EDM4hep

[BrieucF]

Check duplicate issues.

• Checked for duplicates

Goal

In some cases, all tracker hits are attached to the primary particle while there was clearly a secondary produced (colors reflect the MCParticle attached to the trackerHit):
.

This behavior gets cured when enabling --part.keepAllParticles:

but the file size becomes too big for practical purposes...

Wanted behavior: if the particle attached to the hit is not kept in the output file, sth.setParticle(mcp); should not be called, or should be called in a way that does not make it point to an existing particle. See: https://github.com/AIDASoft/DD4hep/blob/master/DDG4/edm4hep/Geant4Output2EDM4hep.cpp

This seems related to (but is not a duplicate of) #471 which I guess was opened for lcio output.

Operating System and Version

Alma9

compiler

gcc 11.4.1

ROOT Version

6.30/06

DD4hep Version

master

Reproducer

source /cvmfs/sw-nightlies.hsf.org/key4hep/setup.sh
ddsim --enableGun --gun.distribution uniform --gun.energy "10*GeV" --gun.particle mu- --numberOfEvents 100 --outputFile ddsim_output_edm4hep.root --random.enableEventSeed --random.seed 10 --compactFile $K4GEO/FCCee/CLD/compact/CLD_o4_v05/CLD_o4_v05.xml --gun.thetaMin "89*degree" --gun.thetaMax "89*degree"
python /afs/cern.ch/user/b/brfranco/work/public/k4geo/draw_hits_cld.py ddsim_output_edm4hep.root
display ddsim_output_edm4hep.png

If you want to directly use the output rootfile, the pathological behavior can be seen on event 87 (starting to count at 0).

Additional context

No response

[andresailer]

Is this particle scattering back, or created in the tracker?

[BrieucF]

It is created in the tracker

[BrieucF]

I had a suspicion that the podio python bindings might be causing this, but I see the same behavior in C++

[andresailer]

The quality bit points to the hit coming from a different particle. So using that bit you can exclude the hits, which is probably what you aim for if you don't want the particle stored?

quality bits: [os......]  o: hit from overlay s: hit from secondary not from the MCParticle associated to it

from the slcio output

 [   id   ] |cellId0 |cellId1 |  position (x,y,z)               |   EDep   |   time   |PDG of MCParticle|   (px,  py, pz)   | pathLength  | Quality 
[snip]
 [00000142] |00000641|00000000|(+5.83e+01, -8.65e+00, -1.87e+00)| 3.31e-05 | 2.31e-01 | 00000000000013  |   unknown         |     n/a     |[ s     ]
        id-fields: (system:1,side:0,layer:5,module:0,sensor:0)

 [00000143] |00000641|00000000|(+5.83e+01, -9.89e+00, -1.31e+01)| 1.23e-04 | 2.71e-01 | 00000000000013  |   unknown         |     n/a     |[ s     ]
        id-fields: (system:1,side:0,layer:5,module:0,sensor:0)

 [00000144] |00123521|00000000|(+5.86e+01, -1.08e+01, -3.83e+01)| 1.03e-06 | 8.54e-01 | 00000000000013  |   unknown         |     n/a     |[ s     ]
        id-fields: (system:1,side:0,layer:5,module:15,sensor:0)

 [00000145] |00123521|00000000|(+5.84e+01, -1.13e+01, -4.23e+01)| 2.64e-06 | 7.40e-01 | 00000000000013  |   unknown         |     n/a     |[ s     ]
        id-fields: (system:1,side:0,layer:5,module:15,sensor:0)

I am not sure the PDGCode should be 13 in this case, I was expecting that to be what the actual particle was. And we don't store the momentum at hit, which might help to understand if we hit some threshold for storing particles.

[andresailer]

Changing

## MinimalKineticEnergy to store particles created in the tracking region
SIM.part.minimalKineticEnergy = 0.001 * MeV

Let's you keep more secondaries, including those that make the hit. Adjust to your liking.

I don't think there is anything wrong in the logic

--------------- print out of MCParticle collection --------------- 

  flag:  0x0
  simulator status bits: [sbvtcls]  s: created in simulation b: backscatter v: vertex is not endpoint of parent t: decayed in tracker c: decayed in calorimeter l: has left detector s: stopped o: overlay

[   id   ]index|      PDG |    px,     py,        pz    | px_ep,   py_ep , pz_ep      | energy  |gen|[simstat ]| vertex x,     y   ,   z     | endpoint x,    y  ,   z     |    mass |  charge |            spin             | colorflow | [parents] - [daughters]

[00000132]    0|        13| 9.85e+00,-1.73e+00, 1.75e-01| 6.30e+00, 1.56e-02, 1.94e-01| 1.00e+01| 1 |[     l  ]| 0.00e+00, 0.00e+00, 0.00e+00| 7.00e+03,-3.33e+01, 1.86e+02| 1.06e-01|-1.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [] - [1,2,3,4,5,6,7,8,9,10,11,12,13,>
[00000133]    1|        11| 9.70e-06, 4.39e-05,-2.79e-06|-0.00e+00, 0.00e+00,-0.00e+00| 5.13e-04| 0 |[s vt  s ]| 2.41e+01,-4.21e+00, 4.28e-01| 2.41e+01,-4.15e+00, 4.24e-01| 5.11e-04|-1.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [0] - [] 
[00000134]    2|        11| 5.23e-04,-4.27e-04,-6.82e-04|-0.00e+00, 0.00e+00, 0.00e+00| 1.09e-03| 0 |[s vt  s ]| 5.86e+01,-1.02e+01, 1.04e+00| 5.89e+01,-1.06e+01,-3.82e+01| 5.11e-04|-1.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [0] - [28,29,30,31,32] 
[00000135]    3|        11|-5.92e-06,-5.05e-05, 7.32e-06| 0.00e+00,-0.00e+00, 0.00e+00| 5.14e-04| 0 |[s vt  s ]| 1.48e+02,-2.53e+01, 2.62e+00| 1.48e+02,-2.54e+01, 2.64e+00| 5.11e-04|-1.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [0] - [] 
[00000136]    4|        11| 9.82e-07,-1.12e-05,-4.32e-05| 0.00e+00,-0.00e+00,-0.00e+00| 5.13e-04| 0 |[s vt  s ]| 3.18e+02,-5.27e+01, 5.62e+00| 3.18e+02,-5.27e+01, 5.56e+00| 5.11e-04|-1.00e+00| 0.00e+00, 0.00e+00, 0.00e+00|  (0, 0)   | [0] - [] 

 [   id   ] |cellId0 |cellId1 |  position (x,y,z)               |   EDep   |   time   |PDG of MCParticle|   (px,  py, pz)   | pathLength  | Quality 
------------|--------|--------|---------------------------------|----------|----------|-----------------|-------------------|-------------|---------
 [00000168] |00122881|00000000|(+1.34e+01, -2.34e+00, +2.37e-01)| 1.32e-05 | 4.54e-02 | 00000000000013  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,side:0,layer:0,module:15,sensor:0)

 [00000169] |00123009|00000000|(+1.44e+01, -2.53e+00, +2.56e-01)| 1.54e-05 | 4.90e-02 | 00000000000013  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,side:0,layer:1,module:15,sensor:0)

 [00000170] |00000257|00000000|(+3.53e+01, -6.15e+00, +6.27e-01)| 1.68e-05 | 1.19e-01 | 00000000000013  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,side:0,layer:2,module:0,sensor:0)

 [00000171] |00000385|00000000|(+3.63e+01, -6.33e+00, +6.46e-01)| 1.33e-05 | 1.23e-01 | 00000000000013  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,side:0,layer:3,module:0,sensor:0)

 [00000172] |00000513|00000000|(+5.73e+01, -9.94e+00, +1.02e+00)| 1.19e-05 | 1.94e-01 | 00000000000013  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,side:0,layer:4,module:0,sensor:0)

 [00000173] |00000641|00000000|(+5.83e+01, -1.01e+01, +1.04e+00)| 2.47e-05 | 1.98e-01 | 00000000000013  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,side:0,layer:5,module:0,sensor:0)

 [00000174] |00000641|00000000|(+5.83e+01, -8.65e+00, -1.87e+00)| 3.31e-05 | 2.31e-01 | 00000000000011  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,side:0,layer:5,module:0,sensor:0)

 [00000175] |00000641|00000000|(+5.83e+01, -9.89e+00, -1.31e+01)| 1.23e-04 | 2.71e-01 | 00000000000011  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,side:0,layer:5,module:0,sensor:0)

 [00000176] |00123521|00000000|(+5.86e+01, -1.08e+01, -3.83e+01)| 1.03e-06 | 8.54e-01 | 00000000000011  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,side:0,layer:5,module:15,sensor:0)

 [00000177] |00123521|00000000|(+5.84e+01, -1.13e+01, -4.23e+01)| 2.64e-06 | 7.40e-01 | 00000000000011  |   unknown         |     n/a     |[   0   ]
        id-fields: (system:1,side:0,layer:5,module:15,sensor:0)

------------|--------|--------|---------------------------------|----------|----------|-----------------|-------------------|-------------|---------

[BrieucF]

Hi Andre, thank you very much for looking into this!

I am not sure the PDGCode should be 13 in this case, I was expecting that to be what the actual particle was

If we decide that the logic is "if the secondary is not kept, the relation points to the youngest kept mother", then PDG could be 13 to be consistent. But this only for lcio, right (edm4hep::SimTrackerHit does not have this member)?

I am indeed ok with the current behavior. I got mislead by the fact that isProducedBySecondary is set to False when the secondary is kept (though the hit does come from a secondary). Depending on how much flexibility one wants to leave to the user of the quality member of edm4hep::SimTrackerHit we may document this behavior here, or in DD4hep.

[andresailer]

though the hit does come from a secondary

No, the hit comes from the particle that caused the hit, not from something it produced. So this sounds technically correct. Otherwise we would essentially always set that bit? If you find the explanation in EDM4hep confusing, you can propose a change there.

[BrieucF]

Yeah, agreed. Thanks Andre!