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( ***************************************************************************
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.
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.
