MONTE CARLO COMPARISONS By: Cristobal Padilla and Matteo Cavalli-Sforza Los Alamos, August 28 1995 SUMMARY: We describe some air shower simulation work we did at LAMPF in August 1995. We found a problem with a version of EAShower that has been widely used within MILAGRO over the last one or two years. This motivated some systematic checks using EGS4 and GEANT 3.21; the energy arriving at MILAGRO (after 750 g/cm**2 of air) from 100 GeV and 300 GeV photon showers was found to be essentially the same when using these two simulation packages, as are the photon and electron multiplicities per shower. These results are in disagreement with analogous simulations using EAShower. Some comparisons are also made with a version of EAShower provided by Tom Yang. On a separate but related issue, we also found that are using EAShower improperly when we simulate muons and heavier particles down to energies comparable to their rest mass. 1. THE SIMULATION PACKAGES (a) EAShower: the first version we used is routinely used at LAMPF, and is identical (but for one format statement) to the version provided to us by Tony Shoup. It can be found in dsirae in the directory /home/padilla/eashower. We also used a version given to us by Tom Yang; however we were later advised by Tom that he has a more recent version, that we did not have the time to run. The one that we used can be found in milagro computer in directory /disk03/people/padilla/easty. In both these packages we chose the energy value ("EGStail" from here on) below which the shower development is simulated with EGS4. (b) EGS4: following a suggestion from D. Williams, we extracted EGS4 from MILAGROSIM. The ausgab and howfar subroutines as well as the main program, called airegs.for, can be found in the directory /disk03/people/padilla/egs in milagro. The file FOR012 contains the PEGS file used in the simulation. (c) GEANT 3.21: we used the version available at LAMPF. That is the version released on March 1994. The GEANT code is called eashower_slic_mixt.for (the files chfile.inc and outinf.inc are also needed) and can be found in /disk03/people/padilla/geant (in milagro) or in /home/padilla/geant (in dsirae). The datacards with the pattern *_realatms_mixt.cards are the input files for the simulation of the exponential atmosphere model, and the ones with the pattern *_bl_mixt.cards the ones for the uniform density atmosphere. 2. AIR COMPOSITION AND DENSITY (a) EAShower: Wrotniak's package uses the "U.S. Standard Atmosphere; we did not change it. The EGS part uses a PEGS file with the following composition (in percentage of atoms, NOT of weight!) Nitrogen 77.93% Oxygen 21.15% Argon 0.47% Carbon 0.02% Hydrogen 0.43% ------------------- Total 100.00% which PEGS obtains X0 = 36.62 g/cm**2 (hence 750 g/cm**2 = 20.48 rad. lengths) Tom Yang's package defines 21 layers of air of exponentially decreasing density, composed of 78.84% NItrogen + 21.16% Oxygen, from which PEGS gives X0 = 37.04 g/cm**2, or 20.23 R.L above MILAGRO. This small difference in composition is not entirely negligible. (b) EGS4: in the EGS4 simulations we used the same PEGS file just described. Reason: we did not have a PEGS code to run, so we used what was available! For the same reason, we only ran EGS4 simulations with a single layer of air with uniform density, with thickness = 750/(1.29 E-3) (STP density in g/cm**3) (c) GEANT 3.21: we defined an air composition as in the above PEGS file for consistency's sake; we did one-layer, uniform density simulations, comparing to EGS4, and 21-layer, exponentially varying density simulations, comparing to the two EASHOWER packages. 3. THE SIMULATIONS (a) THE VARIABLES In all simulations electrons and photons were followed down to a total energy of 1 MeV. We typically ran 1000 showers from 100 or 300 GeV photons coming from zenith. We studied: * EGRD, the energy of all particles in a shower that hit the ground after 750 g/cm**2 of air. * EPART, the individual particle energy | * NG, the number of gammas | at ground level * NE, the number of electrons | * R, the mean radius of all particles at ground level. (b) THE PROBLEM WITH EASHOWER This is what triggered all the following simulations. We found that as we INCREASED the EGS tail energy EGR increases! In other words, as we increase the energy range over which EGS rather that EAShower "makes" the shower, the energy that is NOT absorbed in the air increases. See the following table, EGS tail energy, GeV 1 10 100 Mean EGRD, GeV 1.5 1.9 2.8 Mean EPART, MeV 31 28 29 Mean NG 47 63 90 mean NE 2 3 5 The mean energy per particle remains essentially constant as the energy/shower arriving to the ground increases; the number of particles hitting the ground increases as EGRD increases. This is of course a serious bug. (c) EGS4 vs. GEANT 3.21 - UNIFORM-DENSITY ATMOSPHERE, 100 GeV gammas To see where the "truth" may be we ran both EGS4 and GEANT 3.21; we tend to trust EGS4 better than GEANT, but we could not run EGS in an exponentially- layered atmosphere; so this was the first step: EGS4 GEANT 3.21 Mean EGRD, GeV 0.61+-0.04 0.72+-0.05 Mean EPART, MeV 22.4 24.2 Mean NG 24.5 26.5 mean NE 2.8 3.1 We deem the two results to be in good agreement. The errors on EGRD are rms(EGRD)/sqrt(1000); the distributions of EGRD have a huge rms, typically about 2*mean; they are exponentials with tails (let us remind you that for an exponential distribution, rms=mean). AT THIS POINT, RATHER THAN LOOKING FOR WHERE THE BUG IN EASHOWER IS, WE COMPARED TO AN EASHOWER PACKAGE THAT TOM YANG PROVIDED, WHICH TURNED OUT TO GIVE RESULTS INDEPENDENT OF THE VALUE OF THE EGSTAIL CUT. (d) T. Yang's EASHOWER (100 GeV gammas) We give the results obtained with T.Yang's EASHOWER that he left on the LAMPF disks sometime ago. We unfortunatley lost the output of a run with GEANT 3.21 in our exponentially layered atmosphere, that however agreed with the uniform-atmosphere GEANT result. T.Y, T.Y., EGS tail at EGS tail at .LT. 10GeV .LT. 100GeV Mean EGRD, GeV 0.93+-0.05 0.92+-0.06 Mean EPART, MeV 24 25 Mean NG 36 33 mean NE 4 4 mean R, meter 280 283 Comments: EGRD is 200-300 higher than with EGS or GEANT, but we calculated (using the gamma-function parametrization of showers in the Particle Data Book) that we expect 90 MeV more EGRD due to the fact taht Tom has 0.23 fewer R.L.s in his model of the atmosphere. AS far as EGRD goes, Tom's EAShower is ALMOST in agreement with GEANT and EGS. Also note that, as claimed by Tom, the results are independent of EGS tail energy. However, look at the next point! (e) Simulations with 300 GeV gammas We repeated the above simulations (c) and (d) at 300 geV: * UNIFORM ATMOSPHERE, EGS VS GEANT: EGS4 GEANT 3.21 Mean EGRD, GeV 3.5+-0.2 4.0+-0.2 Mean EPART, MeV 25.2 27.7 Mean NG 125 128 mean NE 16 16 mean R, m 157 172 * EXPO. ATMOSPHERE, GEANT VS. TOM YANG"S EAShower: GEANT T.Y, T.Y., EGS tail at EGS tail at .LT. 10GeV .LT. 100GeV Mean EGRD, GeV 3.7+-0.2 5.4+-0.2 5.7+-0.3 Mean EPART, MeV 25.4 26.5 27.4 Mean NG 128 179 185 mean NE 16 23 23 mean R, m 311 240 237 COMMENTS: * GEANT with Exponential atmosphere gives the same results as with a uniform atmosphere (or, we did not screw up the expo. atm model!) * The quasi-agreement of T. Yang's EAShower with EGS/Geant at 100 GeV appears to have been fortuitous. * We were told by Tom (after presenting all of the above in the regular Monday morning Milagro meeting) that he has a better code at UCSC..... 4. DISCUSSION AND CONCLUSIONS ON E.M. SHOWERS In summary, we could not test a "good" EASHOWER, whereas either EGS4 or GEANT 3.21 can be confidently used, in our opinion, to simulate PHOTON showers (protons and nuclei maybe an other story) Other points brought up at the MILAGRO meeting: * What happens if ONLY EAShower is used to simulate an EM shower? Might it agree with EGS4/Geant 3.21? Unfortunately M.C.S. goes back to Europe tomorrow, however Cristobal will try to answer this question in the remaining time (this week). Answer: Running the version of EASHOWER available in Los Alamos (the one that triggered the problem) the following results are found if EGS is not used at all: Photons Photons 100 GeV 300 GeV Mean EGRD, GeV 0.92 4.5 Mean EPART, MeV 33 36 Mean NG 27 117 mean NE 1.5 9 The conclusion is that they the energies are closer but still above the results from EGS and GEANT. * It appears that many simulations done over the last year or so have too much energy hitting the ground, as a result of the bug we found in EAShower! What are the implications for the effective area of Milagro below 1 TeV? We made plots of the energy and average multiplicity hitting the pond (for central incidence) however, as pointed out by Gus, a full Milagrosim job is needed, because below 1 TeV the showers that trigger (even at the 15 PMT level) are in the high-energy tail of EGRD. So we do not know by how much , nor in which direction the effective area estimates must be revised! * C.P. and M.C.S unfortunately no longer have the PS files with the distributions. The copy of the transparencies that were presented will be mailed to T. ShouP, G. Yodh, and several UCSC people, who have been in touch with C.P. and M.C.S. on these issues and may (?) want to see them; additional copies will be left with Gus S. (please send them to anyone interested, Gus!) 5. ADDITIONAL COMMENTS ON MUONS IN EASHOWER, HADRON_INDUCED SHOWERS, ETC. We asked several people whether the energies that come out of EAShower are KINETIC or TOTAL; we were not totally persuaded of the answers, so we read and figured out what Wrotniak did in the DECAY subroutine for the case of the pion decay to mu + neutrino, and found the following: * the decay in the pion cms is properly and very efficiently simulated (we assume but did not check that the same is true for the other decay chains in DECAY) * however the tranformation to the "lab" system is made in the extremely high energy approximation: beta(decaying particle) = 1, for all angles sin(a+b) = a + b, etc. WROTNIAK HIMSELF WARNS IN THE WRITEUP THAT "BELOW 2 GEV" THE SUBROUTINE IS NOT GOOD. * the energies that come out of DECAY.FOR are total, NOT kinetic energies. If one puts thresholds (as it has been done sometimes) NEAR or BELOW the particle (muon, nucleon) rest energies, the resulting energies and angles are garbage. We do not know whether the muon energy spectra used in past simulations are significantly affected by these observations. It should not be a very big job to change DECAY to follow showers down to the energies for which EASHOWER was not planned to be used. ... if you really want to continue using EASHOWER... but that is an other matter, and we better stop here.