Created attachment 607 [details] investigation of heavy ion problem in hadron calorimeter simulations Starting from GEANT4 10.4 we (NA61/SHINE experiment) have observed strange behavior of heavy ions in the lead/scintillator calorimeter structure. See the slides attached for pictures explaining the problem. Heavy ions start to have less energy deposited than before (in GEANT4 10.3 version). Large discrepancy has been observed now comparing MC and experimental data. Larger Z of ion - bigger the effect.
I have added the new attachment concerning test with GEANT4 10.6 and one event with very low energy in calorimeter. It looks like a lot of neutrons with energies of GeV are killed. Probably this will be a key for solution of the problem with heavy ions as well. In addition, the situation is similar (probably?) to the issue #2047: https://bugzilla-geant4.kek.jp/show_bug.cgi?id=2047 - particles with step length = 0 are disappeared (?)
Created attachment 608 [details] neutrons killed with nKiller
Print out from G4NeutronKiller (I have modified the geant4 code) shows that neutrons killed have big global time of track: Example: * G4Track Information: Particle = neutron, Track ID = 3605, Parent ID = 2112 ********************************************************************************************************* Step# X Y Z KineE dEStep StepLeng TrakLeng Volume Process 0 5 cm 5 cm 54.8 cm 9.79 GeV 0 eV 0 fm 0 fm lead initStep G4NeutronKiller::PostStepGetPhysicalInteractionLength aTrack.GetGlobalTime() = 1.15867e+09 aTrack.GetKineticEnergy() = 9790.24 1 5 cm 5 cm 54.8 cm 9.79 GeV 0 eV 0 fm 0 fm lead nKiller Sergey
Update: I have found that at first ion-ion inelastic process the secondary ion (Pb207) can have big global time: Pb208 ion -> ionInelastic -> Yb161 + Pb207 + ... ... G4HadronicProcess: dynParticle->GetDefinition()->GetParticleName() = gamma, time = 55.2508 G4HadronicProcess: dynParticle->GetDefinition()->GetParticleName() = gamma, time = 55.2508 G4HadronicProcess: dynParticle->GetDefinition()->GetParticleName() = e-, time = 55.2508 G4HadronicProcess: dynParticle->GetDefinition()->GetParticleName() = Yb161, time = 55.2508 G4HadronicProcess: dynParticle->GetDefinition()->GetParticleName() = Pb207, time = 1.15867e+09 G4HadronicProcess: dynParticle->GetDefinition()->GetParticleName() = e-, time = 1.15867e+09 G4HadronicProcess: dynParticle->GetDefinition()->GetParticleName() = gamma, time = 1.15867e+09 2104 5 cm 5 cm 52.9 cm 0 eV 38.8 MeV 3.85 um 16.5 m lead ionInelastic ... Then all particles from this Pb207 after ionInelastic process have about the same global time (1.15867e+09). And neutrons are just "nKill-ed" by this time.
Hint: If in file source/processes/hadronic/management/src/G4HadronicProcess.cc I have made a change in function FillResult like this (set time of secondaries just to global time of parent track): .. //G4double time = std::max(aR->GetSecondary(i)->GetTime(), 0.0) + time0; G4double time = time0; .. then calorimeter response becomes normal (see new attachment)
Created attachment 609 [details] simple fix of time for secondaries
I have managed to reproduce the problem, and I can confirm your hint on the cause of the problem, namely neutrons with a large time which are killed by the neutron killer! And I confirm as well that the problem is present in Geant4 10.4, 10.5 and 10.6. The problem is due to an artifact in the treatment of isomers - artifact that is not present if "RadioactiveDecay" is activated (which is not the case for the FTFP_BERT physics list). Before Geant4 version 10.4, we did not transport isomers, so the ones that were created (internally, not as Geant4 tracks) by hadronic interactions were forced to decay promptly. This means that all tracks produced by hadronic interactions were "prompt" and energetic, and the only way for neutrons to reach a global time of 10 microseconds was to have many, many elastic scatterings and therefore loose their initial energy. Therefore, the neutron killer (that kills neutrons with global time above 10 microseconds) kills only low-energy neutrons, and therefore you do not see events with low energy deposition in the PSD calorimeter. Starting from Geant4 10.4 - and therefore affecting also Geant4 10.5 and 10.6 - we enabled the production of isomers. However, we decided to do it "correctly" only when "Radioactive Decay" is enabled (the rationale was to do things properly only for those users who care about isomers, while not impacting the CPU performance for those who do not). In this case, whenever an isomer is produced, it is transported as any other particle, and therefore it looses quickly its energy by ionization and even if it is produced at high energy, it is very likely to decay at low energy, and therefore its decay products receive only a mild Lorentz boost. Because of this, the neutron killer (that kills neutrons with global time above 10 microseconds) should kill in practice only low-energy neutrons, and therefore you should not see events with low energy deposition in the PSD calorimeter. When "Radioactive Decay" is not enabled - as in FTFP_BERT - we follow the above, correct treatment only for the few isomers with lifetime above a high threshold, set (arbitrarily) to 1000 seconds. For all other isomers (the majority, with a lifetime below this threshold) we decided to do the following rough treatment: we create them, with a "correct" lifetime, but we decay them in the same place where they are born, without transporting them. This implies that, if an isomer is produced at high energy, its decay products receive an artificially strong Lorentz boost, because the parent isomer does not loose energy by ionization, as it should do physically. So, in this way, we can have tracks with large global time (please note that the unit of time in Geant4 is nanosecond!) and large kinetic energies, including of course neutrons, and therefore the neutron killer can kill energetic neutrons and produce artificially events with low energy deposition in the PSD calorimeter. For future versions of Geant4, we will keep the good treatment of isomers when "Radioactive Decay" is activated, otherwise we will disable it completely. Therefore, for future versions (> 10.6), you should not have any issue regarding anomalous (wrong) energy depositions in the PSD calorimeter. For the versions G4 10.4, 10.5 and 10.6, we suggest one of the following 2 ways to bypass the problem of the bad treatment of isomers: 1. Change the following line: fMaxLifeTime = 1000*CLHEP::second; into: fMaxLifeTime = 1.0*CLHEP::microsecond; in the file: G4DeexPrecoParameters.cc in the directory: geant4/source/processes/hadronic/models/de_excitation/management/src/ 2. Activate "Radioactive Decay" on top of the FTFP_BERT physics list. You can do it easily as following, in your main program: #include "G4RadioactiveDecayPhysics.hh" ... FTFP_BERT *thePL = new FTFP_BERT; thePL->RegisterPhysics( new G4RadioactiveDecayPhysics ); ... The solution 1. requires to rebuild Geant4, whereas 2. doesn't (but you need to rebuild your application). Both solutions have a CPU overhead, small for 1., big for 2. Thank you very much for spotting and reporting this problem!
Hello, there is a possibility do not edit Geant4 code but change the parameter. This may be done in user code when /run/initilise is not yet called and only in the master thread: G4NuclearLevelData::GetInstance()->GetParameters()->SetMaxLifeTime(1*CLHEP::microsecond); Similar line exist s in the G4RadioactiveDecayPhysics.cc VI
For both the next release G4 10.7 (including the beta release 10.7.beta expected by the end of June) and the second patch of 10.6 - 10.6.p02 released last week - the problem has been fixed. We have set as default time limit for isomer production to 1 microsecond in all cases, i.e. isomers with half-lives above this threshold are produced in all cases - whether Radioactive Decay is enabled or not. (Before the fix, if Radioactive Decay was not activated, then only isomers with half-lives above 1000 seconds were produced, causing the problem reported by NA61/SHINE.)
