Problem 1824

Summary:	G4ParticleHPPhotonDist::GetPhotons() needs statistical rejection if repFlag is 0
Product:	Geant4	Reporter:	Artem Zontikov <zontikov.a>
Component:	processes/hadronic/models/neutron_hp	Assignee:	dennis.herbert.wright
Status:	RESOLVED FIXED
Severity:	trivial	CC:	Alberto.Ribon, dennis.herbert.wright
Priority:	P5
Version:	10.2
Hardware:	All
OS:	All

Description Artem Zontikov 2016-02-06 23:47:31 CET

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.

Comment 1 Artem Zontikov 2016-03-06 13:46:20 CET

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;
   }
}

Comment 2 Artem Zontikov 2016-03-06 14:00:21 CET

(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

Comment 3 tkoi 2016-06-29 03:02:29 CEST

Thank you for reporting this problem.
Now we are working on this issue.
I will let you know once we have progress.

Tatsumi

Comment 4 dennis.herbert.wright 2019-04-12 23:11:39 CEST

We have adopted your fix and submitted it for testing.  If all goes well, it will be included in Geant4 10.6

Comment 5 Artem Zontikov 2020-01-31 18:16:07 CET

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.

Comment 6 dennis.herbert.wright 2020-02-13 18:44:05 CET

I'm not sure why this fix didn't get in last time, but it's in now.