Created attachment 88 [details] Plots of Geant4 elastic and total cross section for protons and antiprotons compared to data ALICE published a result titled "Midrapidity antiproton-to-proton ratio in pp collisions at $\sqrt{s} = 0.9$ and $7$~TeV measured by the ALICE experiment" which is available at http://arxiv.org/abs/1006.5432 In lines 161-181 they state the antiproton cross section at low momentum (0.45<p<1.05 GeV/c) in GEANT is wrong by about a factor of two. I looked at the following file which seems to be responsible for determining elastic and inelastic cross section: http://www-geant4.kek.jp/lxr/source/processes/hadronic/cross_sections/src/G4HadronCrossSections.cc This appears to be derived from the G3 GHEISHA code. First I plotted the proton-proton and antiproton-proton cross sections obtained from this code along with the data from the PDG. This is attached. Note that what is plotted is the elastic and total cross section; the inelastic is the difference. Geant4 seems to do OK here. The antiproton cross section may be 20% high (for both elastic and inelastic) but not a factor of two. The next obvious place to look is how we get from scattering on protons to scattering on nuclei. This section seems to contain a bug which may be responsible for the reported problem. Line 1426 is where the cross section calculation starts. 1426-1490 is a linear interpolation of the cross section for a given momentum giving the cross sections for scattering on free protons (sigel and sigin). The next step is to adjust for the atomic mass of the target. In lines 1493 and 1494, the correction factors crin and crel are set to unity. In lines 1501-1560, crin and crel are calculated for protons, neutrons, and pions. This is used to account for knowledge of proton and pion scattering on nuclei (Al, Cu, and Pb). Note that in these calculations crin is calculated as basically (correct inelastic cross section)/(scaled up total cross section). So the correction factor is basically inelastic/total (lines 1528 and 1557). This seems to be OK because in line 1562 the correction is applied as crin*total so you get back the inelastic cross section. This is only true, however, for protons, neutrons, and pions which have crin calculated in lines 1528 and 1557. The remaining particles will have an inelastic cross section approximately equal to the total cross section. This affects antiprotons, antineutrons, kaons, lambdas, cascades, and omegas. Looking at the plots attached, one can see that reporting the inelastic cross section as the total cross section combined with the 20% mismatch could produce a factor of two problem.
Additional mail from Kevin: Hello all, This is a response to the emails I received from John and Andrea (both of which are reproduced below). I'm sending this to all the people that appeared in either email. I found the slides Andrea pointed me to interesting although it was unclear what the CMS target was and exactly what the y-axis was (signal response). For my purposes I was less interested in the aftermath of an inelastic interaction, just the cross section. The emails mentioned new improvements in the Geant hadronic model (CHIPS) so I decided to try it out using geant4.9.4.b01. I modified the example novices/N02 to create a 1mm thick Be target followed by five 1mm thick Silicon sensitive detectors, each separated by 4cm (imitation of the first five silicon layers in CMS). Then I fired protons and antiprotons at KE=50 MeV and KE=430 MeV, corresponding to p=310 MeV/c and p=1000 MeV/c. I tried the QGSP_BERT and QGSP_BERT_CHIPS models. I used the default Geant and also a fixed version where I replaced crin=1 with crin=sigin/(sigin+sigel) at http://www-geant4.kek.jp/lxr/source/processes/hadronic/cross_sections/src/G4HadronCrossSections.cc#L1494 I then asked what fraction of the initial particles failed to leave a hit in the last chamber as a surrogate for the inelastic collision rate. At KE=50 MeV (p=310 MeV/c) the fraction interacting were: Protons (all models): 1.66%. pbar QGSP_BERT: 7.49% pbar QGSP_BERT fixed: 4.68% pbar QGSP_BERT_CHIPS: 2.71% At KE=430 MeV (p=1000 MeV/c) the fraction interacting were: Protons (all models): 0.72%. pbar QGSP_BERT: 3.26% pbar QGSP_BERT fixed: 1.67% pbar QGSP_BERT_CHIPS: 1.64% As stated by Andrea (and confirmed by me), CHIPS uses its own inelastic cross section so the fix in G4HadronCrossSections.cc has no effect on the QGSP_BERT_CHIPS results (at least in this energy range). I also tried Lambdas and verified similar behavior. The results from either fix (QGSP_BERT with fix or QGSP_BERT_CHIPS) both show the expected behavior (lower inelastic cross section for pbar). The effect is clearly not small (in this energy range). It would be nice to have the bug fixed for the next release so we have a choice of two reasonable options rather than being forced to use a particular physics model choice. While I have your attention, I wonder if the elastic scaling is also incorrect in the G4HadronCrossSections.cc code. The relevant lines are: 1493 crel = 1.; // this is the default For protons and pions, crel is calculated using interpolation of results on various nuclei like: 1525 crel = xsecel/(0.36*sigel*std::pow(G4double(csa[i]), 1.17)); For all particles, sigel is corrected as: 1563 sigel = crel*0.36*sigel*std::pow(a, 1.17); So for protons and pions, a factor of 0.36 cancels out while for other particles this factor does not cancel out. If one were to actually run this for a=1 (which we don't) one would get an elastic cross section that is too low by 0.36. But since I don't know where this factor comes from, I'm not sure it is a bug. It could be something to deal with differences between protons and neutrons for instance. Thanks everyone for responding, Kevin
Hi Kevin, Thanks very much for your detailed study. It certainly points out a shortcoming in G4HadronCrossSections which we need to address. However, it does not appear to be a bug. As you surmised, this code was derived from the GEISHA code, which for anti-nucleons and kaons, is rather crude due to lack of data to parameterize. The proton-nucleus and antiproton-nucleus cross sections are handled intrinsically differently, although the code makes it look similar due the way it uses the factors "crin" and "crel". At the end of the code, we get, for protons (and neutrons): sigma(inelastic) = sigma(proton_nucleus_inelastic)*(A0/A)**0.685 where sigma(proton_nucleus_inelastic) is the p-nuclear cross section for H, Al, Cu, or Pb, A0 is A for each of these, and A is for the nucleus in question. For antiprotons, we get sigma(inelastic0 = 1.3*sigma(antiproton_proton total)*A**0.63 which is the normal, though too simple, way to go from the elementary reaction to the nuclear one. The factor "crin" therefore changes its meaning from one case to the other in order to get the final forms to come out right. If you ask me, this is confusing programming, but perhaps it has to do with the fact that it started out as Fortan. Anyway, to change "crin" to sigin/(sigin+sigel) for antiprotons would be incorrect, and I'm certain it was not the original author's intent. So, to go forward, some development is needed. The ALICE result gives us a starting point at low energy (Geant4 is too high) while CMS tells us that at high energy Geant4 is too low. If we take a cue from the Glauber-Gribov and CHIPS models, G4HadronCrossSections can be improved by applying an energy-dependent factor to its cross sections, in addition to the usual A scaling, but this will need some study and validation.
