Modeling of rigidity dependent CORSIKA simulations for GRAPES-3

A.Jain; B.Hariharan; B.S.Rao; P.Jagadeesan; P.K.Mohanty; S.Kawakami; S.K.Gupta; S.R.Dugad; S.S.R.Inbanathan; Y.Hayashi

REVIEW 1 major objections 5 minor 24 references

Modeling of rigidity dependent CORSIKA simulations for GRAPES-3

T0 review · 1 major / 5 minor · reviewed 2026-08-14 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read By rejecting primaries below the direction-dependent geomagnetic cutoff, a modified CORSIKA produces the same muon physics in about a third of the simulation time, with half the output data size.

desk verdict Useful CORSIKA Rc pre-filter with a real adoption story, but Table 1's headline speedup numbers are inconsistent with the paper's own energy-scaling assumptions. read the letter →

arxiv 1908.05948 v1 pith:CJUSRZUK submitted 2019-08-16 astro-ph.IM physics.comp-ph

B.Hariharan , S.R.Dugad , S.K.Gupta , Y.Hayashi , S.S.R.Inbanathan , P.Jagadeesan , A.Jain , S.Kawakami

show 2 more authors

P.K.Mohanty B.S.Rao

This is my paper · ORCID

classification astro-ph.IMphysics.comp-ph

keywords cosmicraysgeomagneticcutoffrigidityCORSIKAMonteCarlosimulationairshowerGRAPES-3muontelescopefilter

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

This paper targets a practical bottleneck in air-shower Monte Carlo simulation: for GRAPES-3, more than 60% of randomly generated primary cosmic rays have rigidity below the geomagnetic cutoff of their arrival direction, so they cannot reach the atmosphere, yet the standard CORSIKA simulation tracks them anyway. The authors add a direction-dependent rigidity filter to CORSIKA so that, before shower development begins, each primary's rigidity is compared with a precomputed cutoff-rigidity map and sub-threshold primaries are discarded. In a $10^8$-proton test, this reduced simulation time by a factor of about 3 and output data size by a factor of about 2, while the muon energy spectrum remained identical to that of unmodified CORSIKA. The developers of CORSIKA adopted the feature, and it is present in official versions from v75600 onward. If the claim holds, experiments at locations with large geomagnetic cutoffs can obtain the same physics content with roughly one-third the computing cost.

What carries the argument

The engine of the method is a precomputed geomagnetic cutoff-rigidity map $R_c(\theta,\phi)$, built once per observatory by back-tracing an antiproton through the IGRF-11 geomagnetic field on a 1$^\circ \times$ 1$^\circ$ grid in zenith and azimuth. At initialization, CORSIKA loads this map; for each sampled primary, the cutoff for its direction is obtained by linear interpolation, and if the primary's rigidity is below $R_c$ (within a $-10\%$ tolerance) the event is discarded before interaction tracking. The map is converted from the back-tracing azimuth convention, which runs clockwise from north, to CORSIKA's counter-clockwise convention by a 180$^\circ$ rotation. The filter is implemented for proton and helium primaries, which dominate the composition, and can be extended to other primaries.

What would settle it

Re-run the case A / case B comparison with identical random seeds and extend the comparison beyond muons to the gamma-ray, electron, and hadron spectra at the observational level; a statistically significant deviation in any of these, beyond the stated random-number perturbations, would show that the filter biases the simulated physics. Alternatively, recompute $R_c$ for the GRAPES-3 site with a different geomagnetic field model, such as IGRF-12, and check whether the 56% rejection fraction and the muon spectrum remain materially the same.

Watch

Extended reading notes

Core claim

The central claim is that a precomputed, direction-resolved geomagnetic cutoff rigidity can be used inside CORSIKA to reject forbidden primaries before any shower tracking, with no change to the physics of the accepted events. The paper demonstrates this by generating $10^8$ proton primaries in GRAPES-3's angular and energy range with unmodified CORSIKA (case A) and with the modified code (case B). Case B rejects 56% of primaries before simulation, reduces wall-clock time from 156 to 46 minutes and output file size from 239 to 114 GB, and produces a muon energy spectrum that is statistically identical to case A. Small perturbations at high energies are attributed to the changed random-number sequence caused by skipping rejected events. The paper also argues that the implementation's $-10\%$ tolerance on cutoff rigidity makes the same simulated data reusable for studies in which $R_c$ changes, such as geomagnetic storms.

Load-bearing premise

The load-bearing assumption is that the precomputed cutoff-rigidity map is accurate for every direction the experiment observes; if the map is systematically wrong, the filter will either discard primaries that actually reach the atmosphere or admit primaries that do not, and the paper's muon-spectrum comparison would not reveal a bias localized in the directions or rigidities where the map errs.

Editorial extensions

If this is right

For GRAPES-3, a fixed computing campaign can cover roughly three times as many accepted showers, and storage requirements drop by about half; the paper notes that a thunderstorm simulation that took two months would have exceeded six months and roughly 30 TB of storage without the filter.
Keeping a small tolerance in the rigidity check leaves borderline primaries in the sample, so the same data set can be reused for later analyses where $R_c$ is lower, such as during geomagnetic storms.
Other ground-based cosmic-ray observatories can adopt the method by generating their own back-traced $R_c$ database; the implementation is already part of official CORSIKA from v75600 onward.
Because the filter changes only which primaries enter the shower simulation and not the shower physics itself, all comparisons made on accepted events, such as spectral slopes or directional anisotropies, remain valid.
With zero tolerance, the rejection fraction would rise to about 63%, giving even larger savings at the cost of re-simulation if $R_c$ changes.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

The paper validates the filter with the muon spectrum only; by the logic of the method, gamma-ray, electron, and hadron spectra should also be unaffected, but that is an extension the paper does not demonstrate.
The same prefiltering idea could be applied to other shower codes, such as AIRES or CRY, or wrapped into a common pipeline, wherever a cutoff map is available.
Because the filter reorders the random-number sequence, matching random seeds between filtered and unfiltered runs would isolate whether the small high-energy spectral perturbations are purely numerical or carry any physical content.
The $-10\%$ tolerance creates a tunable trade-off between computational savings and robustness: experiments that expect frequent geomagnetic disturbances should keep a wide tolerance, while those with stable cutoffs can tighten it and save more.

Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

Colleague, here's the short version: the Rc pre-filter is a genuinely useful modification, and its adoption into official CORSIKA v75600 is a real endorsement. But the paper's headline efficiency numbers do not hold up. Table 1 says 156 minutes for 10^8 unfiltered proton showers and 46 minutes for 0.44x10^8 accepted showers in the filtered run. With a -2.7 power law over 10^10-10^13 eV and the stated linear scaling of simulation time with primary energy, the accepted subset is higher-energy, so the expected total time is something like 0.6-0.7 of the unfiltered run, or roughly 94-109 minutes. Even the trivial constant-per-event bound gives 0.44*156=69 minutes. 46 minutes is below that bound. The data size ratio (114/239=0.48) has the same problem: if file size scales with energy, you'd expect ~0.7. So the claimed 3x time reduction and 2x data reduction are not supported by the paper's own assumptions. The timing comparison must be apples-to-apples, but the paper doesn't tell us whether these are CPU or wall-clock times, job sizes, or cluster conditions. This is a central flaw in the main claim.

What the paper does well: the method itself is sound. The Rc map is computed externally by back-tracing with IGRF-11, so no circularity. The comparison of muon energy spectra from filtered and unfiltered sets is appropriate and the spectra agree, which is the right sanity check. The -10% tolerance discussion is a nice practical detail, and the authors honestly cite their ICRC 2015 report, so the novelty is incremental, not invented. The implementation details (coordinate rotation, interpolation) will be useful to anyone who wants to use or adapt the feature.

The soft spots beyond the timing: validation is thin. Only the all-direction muon energy spectrum is compared. I'd want to see directional distributions, muon multiplicity, and maybe shower size before trusting the filter for physics analyses. The Rc map accuracy is not independently checked; the spectrum agreement is suggestive but not a strong test. The high-energy perturbations are waved off as random-number effects, which is plausible but unquantified.

Who is this for? Groups running large CORSIKA production sets, especially at low-energy, geomagnetic-latitude experiments. It's an engineering report, not a physics breakthrough. It deserves a serious referee, because the method is useful and the official adoption is evidence that it works. But Table 1 needs to be corrected or explained. I'd send it to review with a request to fix the timing comparison and to validate the Rc map more directly.

Referee Report

1 major / 5 minor

Summary. The paper presents a modification to the CORSIKA air-shower simulation code that rejects primary cosmic rays whose rigidity is below the direction-dependent geomagnetic cutoff rigidity Rc before shower development, using a precomputed back-traced Rc map for the GRAPES-3 site. The authors report an approximately threefold reduction in simulation time and a twofold reduction in output data size, and they show that the muon energy spectrum from the filtered simulation matches that of the unmodified simulation. The modification has been incorporated into official CORSIKA releases since v75600.

Significance. The work addresses a practical bottleneck for experiments that require billions of simulated primaries: it is a simple, adoptable idea, and its incorporation into official CORSIKA is a strong indicator of usefulness. The validation strategy, comparing final physics output rather than only counting rejected events, is appropriate. However, the quantitative timing claim needs scrutiny, as the reported speedup is not consistent with the paper's own energy-scaling statement.

major comments (1)

[§5, Table 1] The reported simulation times are internally inconsistent with the linear energy scaling stated in §1. For a 10^10–10^13 eV proton sample with an E^-2.7 spectrum, the 44% of primaries that pass the rigidity cut are statistically biased toward higher energies. Since the unmodified run took 156 min for 10^8 showers, a constant-per-event lower bound for the 0.44×10^8 accepted showers is 0.44×156 ≈ 69 min, and any monotone-increasing per-event cost would push the expected time higher (an energy-proportional cost gives roughly 94–110 min). The reported 46 min is below the constant-per-event bound by about 33%, while the data-size comparison (114 GB versus 0.44×239 ≈ 105 GB) is consistent with the accepted subset being higher-energy. This contradiction suggests the two runs are not directly comparable, or that Table 1 contains an error. Because the factor ~3 reduction is a headline claim, please clarify how the times were measured and either correct the table or repeat the benchmark under identical conditions.

minor comments (5)

[Abstract and §5] The abstract states that over 60% of simulated events do not reach the atmosphere, while §5 reports a 56% rejection fraction with the -10% tolerance and estimates about 63% at zero tolerance; these numbers should be reconciled in the text.
[§4] The sentence describing the azimuth conversion from the back-tracking convention to the CORSIKA convention says a 180° rotation is required; converting a clockwise-from-north azimuth to a counter-clockwise-from-north azimuth is a reflection (φ_CORSIKA = 360° − φ_astro), not a 180° rotation, and the text should be corrected or clarified.
[§5, Fig. 4] The claim that the muon energy spectra are identical is supported only by visual inspection; a quantitative comparison, such as a bin-by-bin chi-square or a Kolmogorov-Smirnov test, would make the validation more compelling.
[Table 1] The paper does not state whether the reported simulation times are wall-clock or CPU times, which CORSIKA version was used for the benchmark, or whether the 1000 jobs were run in a single batch; these details should be provided so the efficiency comparison can be reproduced.
[§3] The Rc map is computed with IGRF-11 coefficients but no epoch is specified; because the geomagnetic field changes secularly, the validity period of the map should be stated.

Circularity Check

0 steps flagged · score 0.0 of 10

No circularity found: the Rc cutoff map is an independent back-tracing input, and the claimed speedup and data-size reduction are measured outcomes, not fitted predictions.

full rationale

The paper's central claim is that pre-filtering CORSIKA primaries by direction-dependent geomagnetic cutoff rigidity reduces simulation time and output size without changing the physics content. The Rc database is computed independently by back-tracing antiprotons using IGRF-11 coefficients (Section 3), not from CORSIKA output or from the quantities being compared. The modified code simply rejects a primary before shower simulation when its rigidity is below the interpolated Rc for its arrival direction (Section 4). The claimed factors of roughly 3 and 2 in Table 1 are measured wall-clock times and file sizes, not values fitted to reproduce a target. The physical validation compares muon energy spectra from the modified and unmodified codes; the unmodified case applies the same Rc rejection in post-simulation analysis, so matching spectra is a consistency check of the implementation rather than a tautology. No equation in the paper reduces to an input by construction, and no fitted parameter is renamed as a prediction. References to the authors' earlier work appear only as background physics examples (thunderstorm simulations, geomagnetic-storm muon bursts) and as the original conference report of the method; none is invoked to prove the efficiency claim. A reader could question whether the 46-minute timing is apples-to-apples given the paper's statement that simulation time increases linearly with primary energy and the reported 44% acceptance fraction, but that is a consistency or reproducibility concern, not circularity. Overall, the derivation chain is self-contained and independent of its inputs.

Assumptions & free parameters 1 free parameters · 2 assumptions · 0 invented entities

No new physical entities are introduced. The method adds a pre-filter driven by an externally computed Rc map; the only hand-set parameter is the -10% tolerance. The central claim depends on the accuracy of that map and on the untested assumption that the muon spectrum fully captures the physics equivalence.

free parameters (1)

Rc validation tolerance = -10%
The authors set the filter to accept primaries with rigidity 10% below the nominal cutoff so that transient geomagnetic disturbances or recalculations do not force re-simulation (Section 5). This tolerance affects the reported rejection fraction (56% instead of 63%) and therefore the claimed time and data savings.

assumptions (2)

domain assumption The back-traced IGRF-11 cutoff rigidity map correctly identifies which primaries reach the atmosphere for each direction.
Section 3 constructs Rc via anti-proton back-tracing on a 1 degree grid; Section 4 uses this map to accept or reject every primary. If the map is biased, the filtered set no longer represents the true arriving flux.
domain assumption Removing rejected primaries before simulation does not change the distribution of the accepted showers.
The validated equivalence in Section 5 assumes that the random number sequence in CORSIKA and the absence of rejected-event outputs do not systematically alter the physics; this is checked only through the muon energy spectrum.

how reviews work

0 comments

Cite this review

Pith. "Pith review of Modeling of rigidity dependent CORSIKA simulations for GRAPES-3." pith.science (2026). https://pith.science/paper/CJUSRZUK

@misc{pith2026190805948,
  author       = {Pith},
  title        = {Pith review of: Modeling of rigidity dependent CORSIKA simulations for GRAPES-3},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/CJUSRZUK}},
  note         = {Machine review of arXiv:1908.05948}
}

read the original abstract

The GRAPES-3 muon telescope located in Ooty, India records 4x10^9 muons daily. These muons are produced by interaction of primary cosmic rays (PCRs) in the atmosphere. The high statistics of muons enables GRAPES-3 to make precise measurement of various sun-induced phenomenon including coronal mass ejections (CME), Forbush decreases, geomagnetic storms (GMS) and atmosphere acceleration during the overhead passage of thunderclouds. However, the understanding and interpretation of observed data requires Monte Carlo (MC) simulation of PCRs and subsequent development of showers in the atmosphere. CORSIKA is a standard MC simulation code widely used for this purpose. However, these simulations are time consuming as large number of interactions and decays need to be taken into account at various stages of shower development from top of the atmosphere down to ground level. Therefore, computing resources become an important consideration particularly when billion of PCRs need to be simulated to match the high statistical accuracy of the data. During the GRAPES-3 simulations, it was observed that over 60% of simulated events don't really reach the Earth's atmosphere. The geomagnetic field (GMF) creates a threshold to PCRs called cutoff rigidity Rc, a direction dependent parameter below which PCRs can't reach the Earth's atmosphere. However, in CORSIKA there is no provision to set a direction dependent threshold. We have devised an efficient method that has taken into account of this Rc dependence. A reduction by a factor ~3 in simulation time and ~2 in output data size was achieved for GRAPES-3 simulations. This has been incorporated in CORSIKA version v75600 onwards. Detailed implementation of this along the potential benefits are discussed in this work.

Figures

Figures reproduced from arXiv: 1908.05948 by the authors.

**Figure 1.** Map of Rc in the field of view (FOV) of G3MT of 2.3 mm. A G3MT module consists of 232 PRCs arranged in 4 layers, with alternate layers placed in mutually orthogonal directions which provides an area of 35 m2 . Two successive layers of PRCs are separated by 15 cm thick concrete. The G3MT permits a two-dimensional reconstruction of muon tracks in two vertical orthogonal planes. The vertical separation of two layers o… view at source ↗

**Figure 2.** CORSIKA simulation flow (Modification) [PITH_FULL_IMAGE:figures/full_fig_p005_2.png] view at source ↗

**Figure 3.** Energy spectrum (eV) of PCR events simulated by (a) unmod [PITH_FULL_IMAGE:figures/full_fig_p007_3.png] view at source ↗

Figures from the paper (1 more)

**Figure 4.** Figure 4: Muon energy spectrum (GeV) obtained from (a) unmodified C [PITH_FULL_IMAGE:figures/full_fig_p008_4.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

24 extracted references · 16 canonical work pages

[1]

Mohanty et al., Astroparticle Physics 79 (2016) 23-3 0

P.K. Mohanty et al., Astroparticle Physics 79 (2016) 23-3 0

work page 2016
[2]

Arunbabu et al., Astroparticle Physics 94 (2017) 22- 28

K.P. Arunbabu et al., Astroparticle Physics 94 (2017) 22- 28

work page 2017
[3]

Hariharan et al., Phys

B. Hariharan et al., Phys. Rev. Lett. 122, 105101 (2019)

2019
[4]

D. Heck, J. Knapp, J.N. Capdevielle, G. Schatz, and T. Thou w, Report FZKA 6019 (1998), Forschungszentrum Karlsruhe; available from http://www-ik.fzk.de/corsika/physics description/corsika phys.html

work page 1998
[5]

Monte Carlo Simulation of Proton-induced Cosmic-ray Cas cades in the Atmosphere, UCRL-TM-229452

work page
[6]

Hariharan et al

AIRES: A system for air shower simulations, arXiv:astro- ph/9911331v1 10 B. Hariharan et al

work page arXiv
[7]

Mohanty et al

P.K. Mohanty et al. Phys. Rev. Lett. 117, 171101 (2016); P.K. Mohanty et al. Phys. Rev. D 97, 082001 (2018)

work page 2016
[8]

Hariharan et al., Proceedings of Science PoS(ICRC2015 )448

B. Hariharan et al., Proceedings of Science PoS(ICRC2015 )448

work page

Show all 24 references

[9]

Gupta et al

S.K. Gupta et al. Nucl. Instrum. Methods A 540, 311 (2005); Y. Hayashi et al. Nucl. Instrum. Methods A 545, 643 (2005)

2005
[10]

Mohanty et al., Proceedings of Science PoS(ICRC201 7)357

P.K. Mohanty et al., Proceedings of Science PoS(ICRC201 7)357
[11]

D. F. Smart and M. A. Shea, Adv. Space Res. 36, 2012 (2005)

2005
[12]

Finlay et al., Geophys

C. Finlay et al., Geophys. J. Int. 183, 1216 (2010)

2010
[13]

Pierog et al., arXiv:1306.0121[hep-ph] (2013)

T. Pierog et al., arXiv:1306.0121[hep-ph] (2013)

2013 arXiv
[14]

Kalmykov, S.S

N.N. Kalmykov, S.S. Ostapchenko, and A.I. Pavlov, Nucl. Phys. B (Proc. Suppl.) 52B (1997) 17

1997
[15]

Ostapchenko, Phys

S.S. Ostapchenko, Phys. Rev. D83 (2011) 014018

2011
[16]

Engel, T.K

R. Engel, T.K. Gaisser, P. Lipari, and T. Stanev, Proc. 26th Int. Cosmic Ray Conf. , Salt Lake City (USA), 1 (1999) 415; E.-J. Ahn, R. Engel, T.K. Gaisser, P. Lipari, and T. Stanev, Phys. Rev. D80 (2009) 094003

1999
[17]

W erner, Phys

K. W erner, Phys. Rep. 232 (1993) 87

1993
[18]

Ranft, Phys

J. Ranft, Phys. Rev. D51 (1995) 64; arXiv: hep-ph/9911213 and hep-ph/9911232 (1999)

1995 arXiv
[19]

Drescher, M

H.J. Drescher, M. Hladik, S. Ostapchenko, T. Pierog, and K. W erner, Phys. Rep. 350 (2001) 93 (preprint hep-ph/0007198 (2000))

2001 arXiv
[20]

Fesefeldt, Report PITHA-85/02 (1985),R WTH Aachen; available from: http://cds.cern.ch/record/162911/ﬁles/CM-P00055931.pdf

H. Fesefeldt, Report PITHA-85/02 (1985),R WTH Aachen; available from: http://cds.cern.ch/record/162911/ﬁles/CM-P00055931.pdf

1985
[21]

Fass` o, A

A. Fass` o, A. Ferrari, S. Roesler, P.R. Sala, G. Battisto ni, F. Cerutti, E. Gadioli, M.V. Garzelli, F. Ballarini, O.Ottolenghi, A. Empl and J. Ranft, The physics models of FLUKA: status and recent developments , Computing in High Energy and Nuclear Physics 2003 Conference (C...

2003 arXiv
[22]

Bass et al., Prog

S.A. Bass et al., Prog. Part. Nucl. Phys. 41 (1998) 225; M. Bleicher et al., J. Phys. G: Nucl. Part. Phys. 25 (1999) 1859; http://urqmd.org/

1998
[23]

https://obamawhitehouse.archives.gov/the-press-o ﬃce/2016/10/13/executive-order- coordinating-eﬀorts-prepare-nation-space-weather-events

2016
[24]

Severe Space W eather EventsUnderstanding Societal and Economic Impacts: A W ork- shop Report (National Academies Press, W ashington, 2008) d oi:10.17226/12507; http://www.oecd.org/gov/risk/46891645.pdf

2008

Pith tools

Reviewed August 14, 2026 · model on record in the stance chip above.

[1] [1]

Mohanty et al., Astroparticle Physics 79 (2016) 23-3 0

P.K. Mohanty et al., Astroparticle Physics 79 (2016) 23-3 0

work page 2016

[2] [2]

Arunbabu et al., Astroparticle Physics 94 (2017) 22- 28

K.P. Arunbabu et al., Astroparticle Physics 94 (2017) 22- 28

work page 2017

[3] [3]

Hariharan et al., Phys

B. Hariharan et al., Phys. Rev. Lett. 122, 105101 (2019)

2019

[4] [4]

D. Heck, J. Knapp, J.N. Capdevielle, G. Schatz, and T. Thou w, Report FZKA 6019 (1998), Forschungszentrum Karlsruhe; available from http://www-ik.fzk.de/corsika/physics description/corsika phys.html

work page 1998

[5] [5]

Monte Carlo Simulation of Proton-induced Cosmic-ray Cas cades in the Atmosphere, UCRL-TM-229452

work page

[6] [6]

Hariharan et al

AIRES: A system for air shower simulations, arXiv:astro- ph/9911331v1 10 B. Hariharan et al

work page arXiv

[7] [7]

Mohanty et al

P.K. Mohanty et al. Phys. Rev. Lett. 117, 171101 (2016); P.K. Mohanty et al. Phys. Rev. D 97, 082001 (2018)

work page 2016

[8] [8]

Hariharan et al., Proceedings of Science PoS(ICRC2015 )448

B. Hariharan et al., Proceedings of Science PoS(ICRC2015 )448

work page