Hello. The for-loop inside lines #708-719 selects a photon from the list of discrete photons according to their production cross section (i.e. MF==13) but it misses the reaction cross section itself. After this loop there should be a statistical sampling which returns zero thePhotons pointer if thePartialXsec[iphoton].GetXsec(anEnergy)/theTotalXsec.GetXsec(anEnergy) < G4UniformRand(). In fact, current implementation ignores production cross section, for example, making every O-16(n,na)C-12 reaction coming with 4.438-MeV photon.
Fix: 1. In G4ParticleHPInelasticBaseFS.cc Rewrite line 161 as: theFinalStatePhotons->InitPartials(theData, theXsection); 2. In G4ParticleHPInelasticCompFS.cc Rewrite line 177 as: theFinalStatePhotons[it]->InitPartials(theData, theXsection[50]); 3. In G4ParticleHPPhotonDist.hh Add G4ParticleHPVector* theReactionXsec to private data members, initialize it in the constructor with 0, delete it in the destructor. Rewrite InitPartials as: void InitPartials(std::istream& aDataFile, G4ParticleHPVector* theXsection = 0); 4. G4ParticleHPPhotonDist.cc In InitPartials add the next line: if (theXsec) theReactionXsec = theXsec; After line 719 add the following: if (theReactionXsec) { if (thePartialXsec[iphoton].GetXsec(anEnergy)/theReactionXsec->GetXsec(anEnergy) < G4UniformRand()) { delete thePhotons; thePhotons=0; return thePhotons; } }
(In reply to comment #1) > 3. In G4ParticleHPPhotonDist.hh > Add G4ParticleHPVector* theReactionXsec to private data members, initialize it > in the constructor with 0, delete it in the destructor. > > Rewrite InitPartials as: > void InitPartials(std::istream& aDataFile, G4ParticleHPVector* theXsection = > 0); Correction: there should be theXsec instead of theXsection
Thank you for reporting this problem. Now we are working on this issue. I will let you know once we have progress. Tatsumi
We have adopted your fix and submitted it for testing. If all goes well, it will be included in Geant4 10.6
For some reason version 10.6 still misses item #1 and item #2 from the fix: theFinalStatePhotons->InitPartials(theData, theXsection); in G4ParticleHPInelasticBaseFS.cc theFinalStatePhotons[it]->InitPartials(theData, theXsection[50]); in G4ParticleHPInelasticCompFS.cc Also, for future use, I think that it would be better to skip delete [] theReactionXsec; in ~G4ParticleHPPhotonDist() as the memory allocation does not take place in this class.
I'm not sure why this fix didn't get in last time, but it's in now.
