dCORSIKA - CORSIKA for AMANDA

The report of the progress done in the latest update (ver. 6.016):

Presentations about CORSIKA

 

 

Updates

It is now possible to use different low energy cutoffs (Elow) for different primaries (proportional to the atomic weight of the primaries). As shown in the plot, for each muon producing primary, 99.9% of muons above Elow,mu=238 GeV are produced by primaries with energies above Elow,pri = 1.58 Elow,mu A, where A is atomic weight of the primary. For He and Fe more detailed plots are available. The lowest recorded values of Elow,pri/(A Elow,mu) for He and Fe are 1.51 and 1.41 respectively for 99.9% ratio. These values are lower than 1.58 chosen for the mixture of primaries. For He most values are above 1.58 with only two below (1.51 and 1.57). Values of ratio for Fe are consistently lower than 1.58. However, the number of muons produced by Fe (and other primaries from Li to Mn) above 238 GeV is 2 orders of magnitude smaller (is only 3.6%) than those produced by H and He. Therefore, having merely 99% ratios above 1.58 (which are 1.80 for He and 1.65 for Fe) should be more than enough to maintain the overall ratio above 1.58 for the mixture of primaries when energy cutoffs of primaries are chosen proportional to their atomic weight (just as when energy cutoffs of different primaries are kept the same, which is the case in the above report).

In the above plots poor statistics at high energies does not allow to determine the 99.9% ratio reliably for energies higher than 20 TeV. Extrapolating from the lower energy region, however, one could state that it never gets below 1.58, since the ratio appears to be constant or even rising (the plot for Fe) with energy.

To enable different low energy cutoff feature specify "SPRIC T" in the INPUTS file. Doing so reduces fluxsum from 0.79296 [per meter2 second sr] to 0.35908 (for low energy cutoff for primaries of 376 GeV), increasing the lifetime of mc files 2.2 times. Increasing the low energy cutoff to the highest quoted in the dCORSIKA update report value of 563 GeV reduces the fluxsum further to 0.17706, i.e. increases the lifetime of mc run another 2 times. It is not recommended to increase the low energy cutoff above 432 GeV, however, if one wants to record 99.9% of primary producing muons above their low energy threshold (with less than 0.1% of muons reaching further than the top of B10 not recorded). For this threshold, fluxsum is 0.28156, i.e. only 28% lower (lifetime per run respectively higher) than for the default value of 376 GeV.

To determine the dimensions of the cylinder around the detector, 5 billion primaries (0.5 billion 10x oversampled) with energies above 860 GeV (muon energy cutoff at 430 GeV) were generated with dCORSIKA for the detector with dimensions 800 m x 400 m (length x radius). The detector lifetime of the sample was 1.5 hours. It corresponded to 465947 triggers produced with mass01v001 background simulation (stdkurt ice model). As the dimensions of the cylinder around the detector were changed (by reducing length and radius), the ratio of the number of the events left outside of the cylinder to the total number of triggers was recorded. This ratio is plotted here vs. radius and length. The blue contours correspond to 10-1 to 10-5 of tracks left outside the cylinder. The green contours correspond to the constant AREASUM = pi2 r(r+l) in km2 sr. One should try to select such cylinder dimensions that the AREASUM is minimal for the smallest possible number of tracks left outside the cylinder that produced the trigger. The plot contains 4 corners, but is symmetrical around l=0 (-r is the same as r), so the actual number of tracks outside of the cylinder with chosen dimensions is 1 to 2 times bigger than the contour value in the plot (if a track misses both top and bottom parts for some fixed radius, it misses the whole detector with that radius; however, the track could be contained between the value of the used radius and the maximum radius of 400 m).

For each fixed fraction AREASUM was minimized. Length and radius which minimize it are plotted vs. the fraction here. The values of length and radius are matched values for each fraction and each part of the detector (top and bottom). Values lt=230 m and r=290 m for the top part and lb=240 m and r=300 m for the bottom part of the detector both correspond to no more than 1.5 10-5 fraction of tracks which missed the cylinder of the detector; recalculating for l=480 m and r=300 m of the whole detector (both top and bottom parts) one gets 1.7 10-5 fraction of tracks missed. It is suggested to use values l=600 and r=300 (which correspond to fraction 1.3 10-5), since they are compatible with currently used settings l=800 and r=400 (dCORSIKA files are the same before randomization; and AREASUM only depend weakly on l, it changes by 15% when replacing 480 m with 600 m). Using these new settings increases the lifetime of each run by 1.8 times at the cost of increasing the fraction of missed tracks by 1 order of magnitude.

These plots show the zenith angle and energy distributions of primaries at the trigger level. The highest value of the zenith angle is 83.8 degrees, at which deviation of the surface of the Earth from flat is only 0.12%, i.e. the flat approximation used before for the mass background production was clearly acceptable. The energy spectrum shows a low energy cutoff at just above 860 GeV which indicates that the previously used value of 860 GeV may need to be lowered. A cutoff at high energies (at just below 1.e7 GeV) was present in the INPUTS file and is now removed.

A new mass01v002 run was performed with the new dCORSIKA files generated for detector dimensions l=800 m x r=400 m, and energy cutoffs for primaries and muons recommended in the dCORSIKA update report (Elow,pri=376 GeV, Elow,mu=273 GeV; the different low energy cutoffs for different primaries were not used) for mam and stdkurt ice models. A total of 20 billion primaries (2 billion 10x oversampled) corresponded to detector lifetime of about 1.5 hours. The plots made with these runs (together with those discussed above produced with the older mass01v001 and previously used dCORSIKA based on version v6.014) are summarized in the following table:

MCentriesareaoptimal dimensionsprimary spectraenergy cut
mass01v001 stdkurt dCORSIKA465947areaoptimal dimensionsprimary spectraenergy cut
mass01v002 stdkurt dCORSIKA487722areaoptimal dimensionsprimary spectraenergy cut
mass01v002 mam dCORSIKA287346areaoptimal dimensionsprimary spectraenergy cut
mass01v002 mam pCORSIKA289059areaoptimal dimensionsprimary spectraenergy cut

Even though pCORSIKA uses a different geometrical construction than cylindrical used by dCORSIKA, it is interesting to see how it compares to dCORSIKA, therefore plots produced with pCORSIKA are included in the above table. For these plots 100 mass01v002 files were used (obtained from /data/disk4/mass01/ama@alizarin, run numbers 60001-60100). Even with azimuth angle randomization, 100x oversampling clearly deteriorates the zenith angle and energy distributions which show significant clustering at high zenith angles and energies. Entries are repeated as many as 69 times. With 10x oversampling (no phi randomization) entries are repeated as many as 6 times and clustering is quite a bit less visible. Since primaries in pCORSIKA are not output with their correct coordinates (none of the two lines describing the primary contain correct xy coordinates), it is impossible to determine shower core location. For the above plots coordinates of the first muon in each event were used instead.

All events with primaries (shower cores) missing cylinder with large dimensions (e.g. length=600 m x diameter=600 m) have at least one muon going much closer to the detector than its shower core. However, even if the coordinates of the highest energy muon are randomized inside the cylinder projection instead of those of the shower core, some events will still miss the smaller cylinder (with the dimensions given above), with the main contribution to the trigger coming from a smaller energy muon which passes much closer to the detector. Here are the area and size plots for the mam dCORSIKA run with only the highest energy muon left in each event. The fraction plotted is the fraction of those (highest-energy) muons passing outside the cylinder with given dimensions.

To determine the lowest energy of the primary which can produce a muon triggering the detector, fraction v(E) of such primaries with energies below some energy E was plotted vs. E. The fraction v(E) appears to have a cutoff at some energy Ecut. To find this energy, the fraction v(E) was fit with the 3-parameter function f(E)=(a(E-Ecut))b in the vicinity of the expected cutoff Ecut. Upper bound of the fitted region was varied between v-1(2 10-5) and v-1(2 10-4) and the stability of the fit was found satisfactory. The average value of the fit parameters in the region of fractions between yellow lines was taken for their final values. The above table contains fits for old dCORSIKA (with energy cut for primaries of 860 GeV), new dCORSIKA (for stdkurt and mam ice models) and pCORSIKA. While the region fit was assumed to have cutoff structure, existence of such cutoff at Ecut (value quoted in the plots) is questionable and cannot be inferred from the plots. The fits change considerably when one increases the fitted region. This means that the fits cannot be applied for fractions much different than those used in the fits. Extrapolating the fits to fractions 10-5, 10-6 and 10-7, we get

ice model10-510-610-7
stdkurt676 GeV622 GeV600 GeV
mam711 GeV625 GeV568 GeV

If different low energy cutoffs of primaries (proportional to their atomic weight A) are used, then one would like to know the low energy cutoff divided by A. It is essentially the same as the constant low energy cutoff values given in the above table (given in the table below). To minimize Elow,mu energies of all muons, producing point showers within MMC's active volume (cylinder with length=800 m x diameter=800 m) were translated into the energies that they would need to have to reach the same depth going straight down. The minimum (10-6 fraction) Elow,mu energy is 287 GeV for both mam and stdkurt ice models. Assuming Elow,mu=285 GeV, energies of all primaries were multiplied by Elow,mu/(A Ecut,mu(x)) to determine Elow,pri. The 10-6 fraction is achieved with 498 GeV and 578 GeV for mam and stdkurt ice models respectively.

MCEcutoff,pri/AElow,muElow,pri/A
mass01v001 stdkurt dCORSIKA579 GeV430 GeV600 GeV
mass01v002 stdkurt dCORSIKA631 GeV287 GeV499 GeV
mass01v002 mam dCORSIKA624 GeV287 GeV578 GeV
mass01v002 mam pCORSIKA824 GeV311 GeV724 GeV

The ratios of Elow,pri/Elow,mu is 1.75 for stdkurt and 2.02 for mam ice models. These numbers are greater than 99.9% ratio suggested in the dCORSIKA report (1.58). Since they differ considerably for different ice models, it is suggested to use the value 1.58 from the dCORSIKA report (which is independent on ice models). Using Elow,mu=285 GeV we get Elow,pri/A=450 GeV. This value is considerably lower than Ecutoff,pri/A=600 GeV (taken at 22 GeV lower than the lowest value of Ecutoff,pri/A in the above table).

A more precise (but dependent on the ice model used) calculation of the low energy threshold of muons presented in a table below considers only muons which produce photons (or point showers which produce photons) which hit detector OMs within the trigger window. The cutoff energies are given for 3 cases: 1: all OMs are considered, 2: only OMs with z-coordinate below 400 m, and 3: below 231.5 m are considered. In the first case muons with energies as low as 242 GeV caused OM hits (second OM on string 17 in the example event), which is just above the threshold used (238 GeV), which was chosen with only OMs below 400 m in mind.

ice modelbelow 231.5 mbelow 400 mall hits
10-5 values
mam348 GeV322 GeV247 GeV
stdkurt335 GeV307 GeV246 GeV
lowest recorded (approx. 10-6) values
mam338 GeV311 GeV242 GeV
stdkurt322 GeV287 GeV239 GeV
Elow,mu, 10-6 values
mam312 GeV305 GeV227 GeV
stdkurt298 GeV279 GeV234 GeV

The cutoff of Elow,mu=273 GeV chosen above is clearly acceptable as long as only OMs below 400 m are taken into account. If string 17 or uppermost OMs on strings 11-13 are to be included in the analysis, this setting is no longer sufficient and has to be decreased to at most 227 GeV.

It is now possible to specify the ratio Elow,pri/Elow,mu independently of the low energy primary cutoff as second parameter to the "LOCUT T 1.58" flag in the INPUTS file. Settings recommended for the mass01 background MC are: Ecutoff,pri/A=600 GeV, Elow,mu=273 GeV, LOCUT(2)=1.58, cylinder length=600 m, radius=300 m, l/d=1. These will result in less than 10-5 of triggered events unaccounted for (about 10-6 with l=800 x r=400). The lifetime of a run with 106 primaries is 2.369 seconds, down from 5.394 seconds for corresponding pCORSIKA run. The difference comes from lower primary energy threshold (600 GeV vs. 800 GeV; about (600/800)1.7=0.61 less lifetime) and bigger geometrical area (2.66 km2 sr vs. 1.94 km2 sr; approximately correspond to 10-5 vs. 5 10-5 triggers outside the generation area; these numbers depend on the ice model; for area similar to pCORSIKA dimensions of the cylinder are l=500 m x r=250 m). Execution speed of dCORSIKA run with 106 primaries (no oversampling; on P3 850 MHz) is about 11200 seconds (execution time of a corresponding pCORSIKA run is 25700 seconds).

Due to angle-dependent low energy cut on primaries (set by "LOCUT T 1.58"), dCORSIKA execution time per second of generated lifetime depends only weakly on the low energy cutoff for primaries (Ecutoff,pri), which can therefore be further lowered below the recommended value of 600 GeV. The table below summarizes execution time and lifetime of dCORSIKA run with 105 primaries.

Ecutoff,pri600 GeV500 GeV432 GeV
execution time1451 sec1270 sec1019 sec
lifetime0.2369 sec0.1722 sec0.1333 sec
ratio (exe/life)612573757644

 

Some more (older) links

Some older analysis pages:

My reports/papers/presentations about dCORSIKA:

 

Optimal geometrical dimensions of the detector