Problem 876 - G4PreCompoundModel treats Pi+ incorrectly
Summary: G4PreCompoundModel treats Pi+ incorrectly
Status: RESOLVED FIXED
Alias: None
Product: Geant4
Classification: Unclassified
Component: processes/hadronic/models (show other problems)
Version: 8.1
Hardware: PC Linux
: P2 normal
Assignee: Gunter.Folger
URL:
Depends on:
Blocks:
 
Reported: 2006-07-07 17:22 CEST by mmarino
Modified: 2008-02-12 12:32 CET (History)
0 users

See Also:


Attachments

Note You need to log in before you can comment on or make changes to this problem.
Description mmarino 2006-07-07 17:22:02 CEST
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.
Comment 1 dennis.herbert.wright 2006-07-07 18:46:59 CEST
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?
Comment 2 mmarino 2006-07-07 19:17:59 CEST
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.
Comment 3 mmarino 2006-07-07 19:30:59 CEST
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,
Comment 4 dennis.herbert.wright 2007-01-09 18:02:59 CET
I'm forwarding this to one of the authors of the Binary cascade model.
Comment 5 Gunter.Folger 2008-02-12 12:32:38 CET
Corrected in release 9.1.

Prior to 9.1 binary cascade was not meant to be used for pions.