| Summary: | G4PreCompoundModel treats Pi+ incorrectly | ||
|---|---|---|---|
| Product: | Geant4 | Reporter: | mmarino |
| Component: | processes/hadronic/models | Assignee: | Gunter.Folger |
| Status: | RESOLVED FIXED | ||
| Severity: | normal | ||
| Priority: | P2 | ||
| Version: | 8.1 | ||
| Hardware: | PC | ||
| OS: | Linux | ||
In your physics list have you assigned the Precompound model directly to a pion, or are you using this model as a "back-end" for a high energy model? I'm pretty sure it's used as a "back-end" for a high energy model...I include the physics list code for clarity: G4TheoFSGenerator *theModel = new G4TheoFSGenerator; G4VIntraNuclearTransportModel* theCascade; theCascade = new G4GeneratorPrecompoundInterface; G4ExcitationHandler* theHandler = new G4ExcitationHandler(); G4PreCompoundModel* thePreEquilib = new G4PreCompoundModel(theHandler); theCascade->SetDeExcitation(thePreEquilib); theModel->SetTransport(theCascade); G4QGSModel< G4QGSParticipants > *theStringModel = new G4QGSModel<G4QGSParticipants>; theModel->SetHighEnergyGenerator(theStringModel); G4QGSMFragmentation *theFragmentation = new G4QGSMFragmentation; G4ExcitedStringDecay *theStringDecay = new G4ExcitedStringDecay(theFragmentation); theStringModel->SetFragmentationModel(theStringDecay); theModel->SetMinEnergy(12.0*GeV); theModel->SetMaxEnergy(100*TeV); ... G4PionPlusInelasticProcess* theInelasticProcess = new G4PionPlusInelasticProcess(); theInelasticProcess->RegisterMe( theModel ); pmanager->AddDiscreteProcess(theInelasticProcess); ... (This is essntially the code from the QGSPPiKBuilder. Also, the Binary cascade is registered to theInelasticProcess above (min 0, max 1.2 GeV). I actually believe the offending code would be due to the Binary cascade usage of the PreCompoundModel since the energies are ~19.2 MeV (well below the set min for the QGS). Let me know if you need any clarification. In fact it is in the PreCompoundModel used within the BinaryCascade. A traceback:
#0 G4PreCompoundModel::ApplyYourself (this=0xf9aae40, thePrimary=@0xbfffae70,
theNucleus=@0xfab58a0) at src/G4PreCompoundModel.cc:88
#1 0x08667636 in G4BinaryCascade::ApplyYourself (this=0xf9df548,
aTrack=@0xbfffae70, aNucleus=@0xfab58a0) at src/G4BinaryCascade.cc:150
#2 0x08a275b4 in G4HadronicInteractionWrapper::ApplyInteraction (
this=0xbfffae6b, thePro=@0xbfffae70, targetNucleus=@0xfab58a0,
I'm forwarding this to one of the authors of the Binary cascade model. Corrected in release 9.1. Prior to 9.1 binary cascade was not meant to be used for pions. |
I am receiving the following error: In /common/majorana/simu/geant4.8.1/source/processes/hadronic/models/util/include/G4Fragment.hh, line 264: ===> G4Fragment::SetNumberOfCharged: Number of charged particles can't be greater than number of particles G4HadronicProcess failed in ApplyYourself call for - Particle energy[GeV] = 0.0192 - Material = Lead - Particle type = pi+ *** G4Exception : 007 issued by : G4HadronicProcess GeneralPostStepDoIt failed. *** Fatal Exception *** core dump *** I have traced this to the fact that the PreCompoundModel doesn't handle pi+ appropriately. The above exception comes from the following code excerpt within G4PreCompoundModel (around line 88): 83 84 // Number of Excited Particles 85 anInitialState.SetNumberOfParticles(1+thePrimary.GetDefinition()->GetBaryonNumber()); 86 87 // Number of Charged Excited Particles 88 anInitialState.SetNumberOfCharged(static_cast<G4int>(thePrimary.GetDefinition()->GetPDGCharge()+.01) + 89 static_cast<G4int>(0.5+G4UniformRand())); 90 I'm afraid the above is slightly difficult to read. anInitialState is a G4Fragment. The call to G4Fragment::SetNumberOfParticles() for a pi+ is done with an argument value of 1 (pi+ has 0 baryon number). However, the function call to G4Fragment::SetNumberOfCharged() can take an argument of 1 or 2, depending upon the value of the G4UniformRand() (0->1). If the argument is 2, then G4Fragment::SetNumberOfCharged() will throw an exception: 258 inline void G4Fragment::SetNumberOfCharged(const G4int value) 259 { 260 if (value <= numberOfParticles) numberOfCharged = value; 261 else 262 { 263 G4String text = "G4Fragment::SetNumberOfCharged: Number of charged particles can't be greater than number of particles"; 264 throw G4HadronicException(__FILE__, __LINE__, text); 265 } 266 } I'm pretty sure I've also seen this error in previous version of G4 and so I don't think it is a new error within the most recent version.