Geant4.10.4.p01 I'm looking at high energy (~2.7TeV / nucleon) ions hitting blocks of carbon for LHC collimation and I need electromagnetic dissociation. In the example I use a Xe 131,54 ion on a block of 20cm long, 0.5m wide square of carbon. I found G4EMDissociation but no documentation on how to use it. Although the classes exist, they're not used anywhere in any other code in Geant4 - ie a physics list I could use. I found a talk on a Geant4 course at SLAC on implementing it and followed that: G4HadronInelasticProcess* inelProcIon = new G4HadronInelasticProcess("IonInelastic", G4GenericIon::GenericIon()); G4EMDissociationCrossSection* crossSectionData = new G4EMDissociationCrossSection(); inelProcIon->AddDataSet(crossSectionData); G4EMDissociation* emdModel = new G4EMDissociation(); inelProcIon->RegisterMe(emdModel); G4ProcessManager* pmanager = G4GenericIon::GenericIon()->GetProcessManager(); pmanager->AddDiscreteProcess(inelProcIon); This works for lower energy (351GeV / Xe 131,54 on carbon) but with a lot of warnings about breaking conservation of energy (~1 in 20 events). -------- WWWW ------- G4Exception-START -------- WWWW ------- *** G4Exception : had012 issued by : G4HadronicProcess:CheckResult() Warning: Bad energy non-conservation detected, will re-sample the interaction Process / Model: IonInelastic / EMDissociation Primary: Ru94 (1000440940), E= 248397, target nucleus (6, 12) E(initial - final) = 11244.2 MeV. *** This is just a warning message. *** -------- WWWW -------- G4Exception-END --------- WWWW ------- Is there a bug here? At the LHC energy, Xe 131 54 will have a total energy of ~351TeV. At this point it crashes horribly: ---> Begin of event: 0 G4EnergyRangeManager:GetHadronicInteraction: counter=1, Ek=2669515.134, Material = carbon, Element = Carbon *0* low=100, high=500000 In /Users/nevay/physics/packages/geant4.10.04.p01/source/processes/hadronic/management/src/G4EnergyRangeManager.cc, line 129: ===> GetHadronicInteraction: No Model found -------- EEEE ------- G4Exception-START -------- EEEE ------- *** G4Exception : had005 issued by : G4HadronicProcess::PostStepDoIt In /Users/nevay/physics/packages/geant4.10.04.p01/source/processes/hadronic/management/src/G4EnergyRangeManager.cc, line 129: ===> GetHadronicInteraction: No Model found Target element Carbon Z= 6 A= 12 Unrecoverable error in the method ChooseHadronicInteraction of IonInelastic TrackID= 92 ParentID= 1 Sb114 Ekin(GeV)= 304325; direction= (-1.10045e-06,1.01819e-06,1) Position(mm)= (-4.72925e-05,4.37574e-05,72.1658); material carbon PhysicalVolume <c1_collimator_pv> No HadronicInteraction found out *** Fatal Exception *** core dump *** -------- EEEE -------- G4Exception-END --------- EEEE ------- *** G4Exception: Aborting execution *** I feel that even if it's beyond the energy range of the model, it should not crash. From the table print out at the beginning, I see: ==================================================================== HADRONIC PROCESSES SUMMARY (verbose level 1) --------------------------------------------------- Hadronic Processes for GenericIon Process: ionInelastic Model: EMDissociation : 100 MeV/n ---> 500 GeV/n Cr_sctns: Electromagnetic dissociation: 0 eV ---> 100 TeV Cr_sctns: GheishaInelastic : 0 eV ---> 100 TeV ================================================================ Is the 100TeV / particle or per nucleon? Am I using the right model and if so am I using it correctly? Any advice would greatly appreciated.
To answer your first question: the warning is just telling you that the model is re-sampling the interaction because it found non-conservation at some level. For a given event, the code will loop until finds an energy-conserving solution. As long as the loop doesn't go on for too long (no problem from what you've shown me) everything is OK. To solve the crash, you need to raise the maximum model energy. By default, it's set to 500 GeV per nucleon because that's our best guess at the upper end of the model's validity. When this is exceeded, the physics list looks for a higher energy model, doesn't find it, and then crashes. This is intentional in order to guarantee that users supply models for all energies likely to occur. In your physics list, after you've instantiated the model, add the line emdModel->SetMaxEnergy(X*TeV); where X is as big as you need it to be. There is no theoretical maximum for this model; it's just that it has not been validated at the TeV scale.
