Problem 178

Summary: Angular and transverse momentum distributions are asymmetric, and related bugs.
Product: Geant4 Reporter: dbailey
Component: processes/hadronicAssignee: Hans-Peter.Wellisch
Status: CLOSED FIXED    
Severity: normal    
Priority: P2    
Version: 2.0   
Hardware: PC   
OS: Linux   

Description dbailey 2000-12-14 05:41:28 CET 
If a hadron travelling along the x axis interacts, the shower
is not symmetric in y. Looking at low energy interactions a
very clear bias can be seen for events in which there is a net
momentum in the +y direction. After some searching this is traced
to G4ReactionDynamics::Rotate, based on the Gehisha fortran rountine defs.F.
It turns out that this bug was found in GEANT3/Gheisha last year by Andrea
Venturi
(see http://www.pi.infn.it/~venturi/Gheisha.html ) and fixed
last November (see /cern/2000/src/geant321/gheisha/defs.F).

     While studying this bug, I also notice the following problems:

***************************************************************************
**** used cos(rthve) when it should be cos(phinve) ****
$g4install/source/processes/hadronic/util/srcg4reactiondynamics.cc
1837:             targetparticle.setmomentum( pp*sin(rthnve)*cos(rthnve)*mev,
1857:             currentparticle.setmomentum( pp*sin(rthnve)*cos(rthnve)*mev,
***************************************************************************

***************************************************************************
**** an isotropic angular distribution is generated by flat
distributions in cos(theta) and phi. this is correctly done many times, but
there are also many places where a uniform theta distribution has been used
instead of a uniform cos(theta) distribution. ****

$g4install/source/processes/hadronic/util/src/g4reactiondynamics.cc
539:                 g4double rthnve = pi*g4uniformrand();
646:       g4double rthnve = pi*g4uniformrand();
749:               g4double rthnve = pi*g4uniformrand();
795:                 g4double rthnve = pi*g4uniformrand();
844:         rthnve = pi*g4uniformrand();
1106:         rthnve = pi*g4uniformrand();
1124:         rthnve = pi*g4uniformrand();
1141:           rthnve = pi*g4uniformrand();
1835:             rthnve = pi*g4uniformrand();
1855:             rthnve = pi*g4uniformrand();
2801:         rthnve = pi*g4uniformrand();
2830:         rthnve = pi*g4uniformrand();
2861:           rthnve = pi*g4uniformrand();
2887:           rthnve = pi*g4uniformrand();
2906:           rthnve = pi*g4uniformrand();
2926:             rthnve = pi*g4uniformrand();

$g4install/source/processes/hadronic/hadronic source
(4.2.0r2):util:src:g4nucleus.cc
263:     g4double theta=randflat::shoot((hepdouble)0.,(hepdouble)pi);  //
isotropic decay

$g4install/source/processes/hadronic/util:src:g4reactionkinematics.cc
49:    g4double theta=randflat::shoot(hepdouble(0.),hepdouble(pi));  // isotropic
decay

$g4install/source/processes/hadronic/
models:neutron_hp:src:g4neutronhpphotondist.cc
272:           g4double theta = pi*g4uniformrand();
370:           g4double theta = pi*g4uniformrand();

$g4install/source/processes/hadronic/stopping:src:g4muminuscapturecascade.cc
108:    g4double theta = (2.0 * g4uniformrand() - 1.0) * pi ;

$g4install/source/processes/hadronic/stopping:src:g4piminusstopmaterial.cc
125:       g4double theta1 = generateangle(pi);
***************************************************************************

***************************************************************************
**** there are several places where 0 < phi <2pi is desired, but the attempt
to fix phi<0.0 values actually pushes them to phi>2pi.  (i.e. since phi is
negative in this case, subtracting it from twopi gives a number bigger than
twopi.) i assume the real intent is to add 2pi to phi. ****

$g4install/source/processes/hadronic/models/high_energy/src/g4heinelastic.cc
1146:                               if( phi < 0.0 )phi = m_2pi - phi;
1175:                               if( phi < 0.0 )phi = m_2pi - phi;
1254:             if( phi < 0.0 )phi = m_2pi - phi;
3379:                               if( phi < 0.0 )phi = m_2pi - phi;
3408:                               if( phi < 0.0 )phi = m_2pi - phi;
3479:             if( phi < 0.0 )phi = m_2pi - phi;

$g4install/source/processes/hadronic/util:src:g4reactiondynamics.cc
461:               if( phi < 0.0 )phi = twopi - phi;
489:               if( phi < 0.0 )phi = twopi - phi;
590:         if( phi < 0.0 )phi = twopi - phi;
719:               if( phi < 0.0 )phi = twopi - phi;
***************************************************************************

***************************************************************************
**** there is no fermi motion for hydrogen, so should subtract 1 from the
finalstatenucleons/targ count in the following lines. ****

$g4install/source/processes/hadronic/util/src/g4reactiondynamics.cc
2685:     ran1 = a1*sin(rx)*0.020*numberoffinalstatenucleons*gev;
2686:     ran2 = a1*cos(rx)*0.020*numberoffinalstatenucleons*gev;

$g4install/source/processes/hadronic/models:high_energy:src:g4heinelastic.cc
1448:    pvmx[6].setmomentum( pvmx[6].getmomentum().x()+ran1*0.020*targ,
1449:                         pvmx[6].getmomentum().y()+ran2*0.020*targ );
2647:    pvmx[7].setmomentum( pvmx[7].getmomentum().x()+ran1*0.020*targ,
2648:                         pvmx[7].getmomentum().y()+ran2*0.020*targ );
3697:    pvmx[6].setmomentum( pvmx[6].getmomentum().x()+ran1*0.020*targ,
3698:                         pvmx[6].getmomentum().y()+ran2*0.020*targ );
4540:    pvmx[7].setmomentum( pvmx[7].getmomentum().x()+ran1*0.020*targ,
4541:                         pvmx[7].getmomentum().y()+ran2*0.020*targ );

***************************************************************************

***************************************************************************
**** the checkmomentum method in g4reactiondynamics seems to require all
final state particle momenta to have a magnitude less than the incoming
particle momentum. this actually causes energy and momentum non-conservation
is cases such as two-body backward scattering of pions off nucleons, since
the forward going nucleon in such scattering has a momentum greater
than the incoming pion momentum. checkmomentum is called after 2 of the 3
invocations of addblacktrackparticles. i don't know why it is not called
the last time. it is not clear to me why checkmomentum is
needed at all, so for the moment i have simply commented it out. ****

$g4install/source/processes/hadronic/util/src/g4reactiondynamics.cc
1201:       momentumcheck( modifiedoriginal, currentparticle, targetparticle,
vec, veclen
2033:       momentumcheck( modifiedoriginal, currentparticle, targetparticle,
vec, veclen
3087:  void g4reactiondynamics::momentumcheck(

***************************************************************************
Comment 1 Hans-Peter.Wellisch 2001-01-23 03:23:59 CET 
Hi,

  thanks for submitting this one. I have looked through it, and it turns out
that most of what you found is actually correct as is, and does not need
changing. The G4ReactionDynamics::Rotate and G4muminuscapturecascade.cc one on
the other hand is valid, and will be fixed.

Many greetings,

Hans-Peter.
Comment 2 Hans-Peter.Wellisch 2002-02-12 12:21:59 CET 
To close this, let me go through the various blocks in some detail:

Block 0: NO PROBLEM Still, this was not meant like this, but there is no harm.
It is contained in an exceptional condition, and invoked only for energies below
1keV.

Block 1: these are indeed lines of code that would lead to phi asymmetries,
if they were not again exception cases. Here the algorithm basically failed to
converge, the secondary energy ending up at a very small value and any direction
will do.

Blocks 2,3: this were bugs indeed. They affects the nucleon distribution of a
not recommended model for K- absorption at rest; was fixed.

Blocks 4,5: important bugs, was fixed at the time the bug report was submitted.

Block 6: Bug affecting a particular (not normally recommended) pi- at rest
model, fixed.

Block 7: This one is with Harm. Since he has not acted on it, I assume it is the
way, he wanted it.

Block 8: Yes, but the error (~1-5MeV) is much smaller than other approximations
made. A rigorous treatment would be much more involved. Here no change.

Block 9: dito.

Block 10 (last): I agree with this, but note that if you are sensitive to the
decription of rare processes like quasi-elastic backward scattering of O(GeV)
pions, this is not the model you should use. Most likely you will want to have
your own custom process for the rare channel of your interest.

Many greetings,

Hans-Peter.