Problem 1680

Created attachment 289 [details]
patch relative to ref-08

This can be fixed by rejecting proposed interactions in G4CascadeInterface. If
- the incident particle is a gamma or electron
- with kinetic energy < 50 MeV
- the outgoing particles contain only gammas
- maximum retries not reached
the proposed interaction is rejected.

This produces the correct number of photoneutrons and should have minimal effect on anything else.

Thanks for your prescription.  We'll give it a try and let you know how it goes.   I have one concern: for light nuclei we seem to be getting good agreement as things now stand.  So we'll check this under the new conditons

I am implementing this special selection as part of G4InuclCollider, just after the nuclear fragment de-excitation (so that it can see both direct cascade secondaries as well as evaporation products).

In order to evaluate performance, it would be very helpful to have some data with which to compare the output of BERT without and with this special selection.  If you have any tabulated data, e.g., for the gamma-W channel, could you provide it (either via e-mail or as an attachment here)?

Daren, you may find a tag of BERT (source/processes/hadronic/models/cascade) in the SVN repository, hadr-casc-V10-01-01, which implements your proposed accept/reject filtering for low-energy photonuclear interactions.  The implementation is *disabled* by default, so that we can evaluate its effect without breaking the reference builds.  

If you wish to test it with your own application and data, set the environment variable G4CASCADE_CHECK_PHOTONUCLEAR = 1 at runtime (not compile time).  I would certainly appreciate any feedback!

  -- Michael Kelsey

It looks to give correct results. Thanks. Some comments that may or may not be relevant:
- I hope it can be implemented by default (without environment variable)
- The 50 MeV limit is arbitrary. It will result in a step in a graph of neutron yield vs. photon energy.
- Hadr03 runs about half the speed with the patch.
- Without the patch, the target nucleus is listed as a generated nucleus. With the patch, it's not. (maybe this simplifies the test?)

Thank you for the followup and summary, Daren.  You wrote:
>> - I hope it can be implemented by default (without environment variable)

Our intention is to do so, once it's validated.  I wanted to be able to put the code into SVN (rather than passing around uncontrolled .cc/.hh files) without affecting any existing runs.

>> - The 50 MeV limit is arbitrary. It will result in a step in a graph of neutron yield vs. photon energy.

This is something which needs to be nailed down.  I took that cut from you're bug report :-)  I think we'll need some guidance from low-energy experts about where that threshold should be, whether it is (A,Z) dependent, etc.  If we can avoid steps, that will always be better.

>> - Hadr03 runs about half the speed with the patch.

Do you mean with the patch enabled (i.e., the envvar), or do you mean even without the gamma-N filter enabled?

>> - Without the patch, the target nucleus is listed as a generated nucleus. With the patch, it's not.
>>    (maybe this simplifies the test?)

That's weird.  The final state should not be any different, provided the cascade passes the filter.

> >> - The 50 MeV limit is arbitrary. It will result in a step in a graph of neutron yield vs. photon energy.
> 
> This is something which needs to be nailed down.  I took that cut from you're
> bug report :-)  I think we'll need some guidance from low-energy experts about
> where that threshold should be, whether it is (A,Z) dependent, etc.  If we can
> avoid steps, that will always be better.

Yes. 50 MeV is above the giant dipole resonance so most of photoneutron production is correct this way. I don't know what happens above that energy. 

> >> - Hadr03 runs about half the speed with the patch.
> 
> Do you mean with the patch enabled (i.e., the envvar), or do you mean even
> without the gamma-N filter enabled?

I was comparing with and without the envvar.

> >> - Without the patch, the target nucleus is listed as a generated nucleus. With the patch, it's not.
> >>    (maybe this simplifies the test?)
> 
> That's weird.  The final state should not be any different, provided the
> cascade passes the filter.

The filter rejects the reaction (e.g.)
 gamma + W186 --> N gamma + W186
All other observed reactions change the target isotope (which makes sense but I'm not an expert).

>> > >> - Hadr03 runs about half the speed with the patch.
>> I was comparing with and without the envvar.

Okay.  Running the BERT internal unit test (which has no geometry, tracking, or process
manager overhead), there is definitely a difference.  For 50k interactions, 10.1 takes 46 seconds.  Turning on the photonuclear check increases that to 66 seconds, a 43% time increase.

>> The filter rejects the reaction (e.g.)
>>  gamma + W186 --> N gamma + W186
>> All other observed reactions change the target isotope (which makes sense but 
>> I'm not an expert).

Right; I misunderstood your comment.  Implementing a comparison of the final nucleus is definitely more efficient than looping over the full final state, though it may not be significant (since the number of final-state gammas is small).  Making this change reduces my test job to 58 seconds, still a 26% increase over not doing the test at all.

I don't know why Hadr03 doubles the execution time.  Maybe it's because there are more hadronic secondaries to be tracked?

Hi Mike, Daren,

have you tried out BERP option when G4 native pre-compound model is used? Or the problem exist on level of the cascade?

I would think that increase of generation CPU is not important if the physics become correct. CPU should be controlled in applications: FullCMS, SimplifiedCalorimeter, medical....

Vladimir

I believe all the physics lists and factories use the same model for photonuclear. I get:

                           Hadronic Processes for gamma

  Process: photonNuclear
        Model:            BertiniCascade: 0 eV  ---> 3.5 GeV


> have you tried out BERP option when G4 native pre-compound model is used? Or
> the problem exist on level of the cascade?

Are you suggesting I modify an existing list to change this line? To what?
  this->RegisterPhysics( new G4EmExtraPhysics(ver) );

I agree the time issue is not relevant.

> have you tried out BERP option when G4 native pre-compound model is used? Or
> the problem exist on level of the cascade?

Using the precompound model returns
===> BAD primary type in G4PreCompoundModel: gamma

Mike's suggestion was to try CHIPS model in Geant4 version 9.5.patch2. It gives realistic results. 

For 12 MeV incident gamma on G4_W, almost every photonuclear interaction should produce 1 neutron.

I modified FTFP_BERT to not have electromagnetic physics and added a UserSteppingAction to Hadr00 example to print the secondaries, step length, and abort the event after first step. For 10,000 incident gamma, I see a step length of 54.5 +- 0.5 cm and 10,000 neutrons produced. 

This is in comparison to 10.1:  MeanFreePath:  54.684 cm  +- 54.272 cm 
crossSection per atom: 289.25 mbarn
   neutrons produced:    2927

The cross section is good for each version. CHIPS in 9.5 gets the channels correct (at least at this energy), Bertini cascade in 10.1 gets the channels incorrect.

This crashes in 10.2, using e.g. Hadr03, gamma.mac, and envvar set.

Here: 
    at /geant4.10.02/source/processes/hadronic/models/cascade/cascade/src/G4InuclCollider.cc:360  

  G4double mfinalNuc  = checkOutput.getOutgoingNuclei()[0].getMass();

The vector returned by checkOutput.getOutgoingNuclei() has no elements.

I did a comparison of some versions of GEANT plus FLUKA for neutron production from muons. Also played with this fix.

Method was taken from arXiv:1302.4275. Neutrons are produced by 260 GeV muon beam going through axis of lead cylinder 3m x 3m x 3m. Calculated is neutron production rate taken as number of born neutrons in central 1.5m section per original muon divided by 1.5m and lead density. Number of neutrons was corrected for neutrons disappearing in NeutronInelastic. Physics list was QGSP_BERT.

Neutron yield per production process in units of 10^-3 neutrons/muon/(g/cm^2)
Process     FLUKA     G4.9.5    G4.10.0   G4.10.0+fix G4.10.2   G4.10.2+fix
gamma-N     1.99614   1.97885   0.78839   1.90095     0.78755   1.86966
neutron-N   1.72939   1.51601   1.59747   1.70263     1.46804   1.60267
pion-N      0.63214   0.46777   0.42166   0.47598     0.41134   0.46071
muon-N      0.25495   0.21250   0.17839   0.17638     0.15345   0.15754
proton-N    0.14301   0.18236   0.17582   0.18778     0.17176   0.19041
other       0.04437   0.16452   0.11842   0.13683     0.12747   0.15140

As one can see the patch actually works and moves photo-nuclear production to a (more or less) right place. Also I plotted an energy distributions of initial energy for born neutrons and of final energy for neutrons escaped from the lead block. Both are attached to this page as (lead-neutrons-yield.png and lead-neutrons-escaped.png). I don't see there 'may be' spectra violations.

Created attachment 377 [details]
Energy spectra of born neutrons

lead-neutrons-yield.png

Created attachment 378 [details]
Energy spectra of escaped neutrons

lead-neutrons-escaped.png

(In reply to comment #12)
> This crashes in 10.2, using e.g. Hadr03, gamma.mac, and envvar set.
> 
> Here: 
>     at
> /geant4.10.02/source/processes/hadronic/models/cascade/cascade/src/G4InuclCollider.cc:360 
> 
>   G4double mfinalNuc  = checkOutput.getOutgoingNuclei()[0].getMass();
> 
> The vector returned by checkOutput.getOutgoingNuclei() has no elements.

I didn't observe crashes during my testing. Though in v10.1 they fixed a bug with improper brackets and that affects exactly the place of the patch so I had to modify it accordingly.

Re #12:  There are conditions under which an interaction can completely disrupt a nucleus.  Obviously, this would count as a true "photonuclear" event, so I've extended the checks in G4InuclCollider::photonuclearOkay() such that if there is no outgoing nucleus, the interaction is good.

I have also removed the envvar "protection," so that photonuclearOkay() is not checked always.  This is driven by V. Belov's recent comments and supporting plots.  This will now be able to go through G4's larger-scale validations, including calorimeter response tests, to ascertain whether it has unexpected side effects.

Gentlemen, 
I am sorry for invading into the discussion.

Lately I did some comparisons between Geant4 10.1.p02 and Fluka 2011.2c.2. In brief, we were simulating neutron source based on e-linac (Ee= ~30MeV). There were two types of target - Tungsten and U-238. Both of them were surrounded with water moderator.

Geant4 and Fluka give the same neutron yield for U-238 while the results for Tungsten were differ by factor of ~2.8. Both codes use practically the same cross-section values for Uranium. Can't say anything about Tungten because I did not managed to extract cross-section from fluka database (it has binary format without description).

As you may know, there are some databases with evaluated data - www.nndc.bnl.gov/sigma . These databases contain cross-section and energy-angle distributions of outgoing particles for photonuclear reactions in ENDF format. Some of them are tabulated up to 140-150 MeV, others only up to 20-30 MeV.

So, it could be possible to combine G4CascadeInterface with a model based on evaluated data with smooth transition in the range of 140 MeV or 20 MeV. The code for handling evaluated data and generating final state could be snipped from particle_hp directory.

Hi! Is this fix already included in version 10.02.p01?
From the release note of 10.02, I understand that it's included, but one needs to set the variable G4CASCADE_CHECK_PHOTONUCLEAR  to 1. Is that the case also in geant4.10.02.p01?
Thanks

Paola

Hi Paola, the fix is present in 10.2.patch 1. The environment variable needs to be set. 

Also there is still a bug. In line 360 of G4InuclCollider.cc

http://www-geant4.kek.jp/lxr/source/processes/hadronic/models/cascade/cascade/src/G4InuclCollider.cc

there needs to be protection against vectors of size 0.

// add this line
  if (checkOutput.getOutgoingNuclei().size() > 0) {
// before this one
    G4double mfinalNuc  = checkOutput.getOutgoingNuclei()[0].getMass();

I'm reopening this bug.

Thanks, Daren.
To fix this in my own local distribution, can you confirm that the code after the fix should be like this:

// Hadron production changes target nucleus
  if (checkOutput.getOutgoingNuclei().size() > 0) {
    G4double mfinalNuc  = checkOutput.getOutgoingNuclei()[0].getMass();  
    G4double mtargetNuc = interCase.getTarget()->getMass();
    if (mfinalNuc != mtargetNuc) return true;	// Mass from G4Ions is fixed
  }
  
  if (verboseLevel>2)
    G4cout << " photonuclear produced only gammas.  Try again." << G4endl;

  return false;		// Final state is entirely de-excitation photons

In particular, are the brackets ok?

Thanks

Paola

Hi Paola, 

Ideally one of the hadronics experts will comment. The problem is occasionally the vector checkOutput.getOutgoingNuclei() will have size zero. Thus we can't just take the [0] element. However in such a case I don't know whether to return true or false. For my work the vector has size zero infrequently and it doesn't matter. 

The code I use is below. Note it will still crash if verboseLevel>2.

G4bool G4InuclCollider::photonuclearOkay(G4CollisionOutput& checkOutput) const {
  if (interCase.twoNuclei()) return true; // A-A is not photonuclear

  G4InuclElementaryParticle* bullet =
    dynamic_cast<G4InuclElementaryParticle*>(interCase.getBullet());
  if (!bullet || !(bullet->isPhoton() || bullet->isElectron())) return true;

  if (verboseLevel>1)
    G4cout << " >>> G4InuclCollider::photonuclearOkay" << G4endl;

  if (bullet->getKineticEnergy() > 0.050) return true;

  if (verboseLevel>2) {
    G4cout << " comparing final nucleus with initial target:\n"
     << checkOutput.getOutgoingNuclei()[0] << G4endl
     << *(interCase.getTarget()) << G4endl;
  }

  // Hadron production changes target nucleus
  if (checkOutput.getOutgoingNuclei().size() > 0) {
    G4double mfinalNuc  = checkOutput.getOutgoingNuclei()[0].getMass();
    G4double mtargetNuc = interCase.getTarget()->getMass();
    if (mfinalNuc != mtargetNuc) return true; // Mass from G4Ions is fixed

    if (verboseLevel>2)
      G4cout << " photonuclear produced only gammas.  Try again." << G4endl;
  }
  else return true;

  return false;   // Final state is entirely de-excitation photons
}

A new photo-nuclear model G4LENDorBERTModel has been developed to fix this problem.  It replaces the Bertini cascade by using photo-nuclear data from the GND database as accessed by the G4LEND models. It covers 167 nuclides (the most common ones) up to 20 MeV.  If a nuclide is encountered which is not in this list, or the energy is above 20 MeV, the model automatically switches back to Bertini.

The model is now being validated and will be available with the Geant4 10.4 release in December.