diff --git a/analysis/do_showers.sh b/analysis/do_showers.sh
new file mode 100755
index 0000000000000000000000000000000000000000..7b9d018139320d4d8a808ce69141930fab4f6825
--- /dev/null
+++ b/analysis/do_showers.sh
@@ -0,0 +1,311 @@
+#!/bin/bash
+# /************************************************************************/
+# /* */
+# /* Package: ARTI */
+# /* Module: do_showers.sh */
+# /* */
+# /************************************************************************/
+# /* Authors: Hernán Asorey */
+# /* e-mail: hernanasorey@cnea.gob.ar */
+# /* */
+# /************************************************************************/
+# /* Comments: Script to automatize the analysis of simulated showers */
+# /* */
+# /************************************************************************/
+# /*
+#
+#
+# Copyright 2021
+# Hernán Asorey
+# ITeDA (CNEA), Argentina
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are
+# met:
+#
+# 1. Redistributions of source code must retain the above copyright
+# notice, this list of conditions and the following disclaimer.
+#
+# 2. Redistributions in binary form must reproduce the above copyright
+# notice, this list of conditions and the following disclaimer in the
+# documentation and/or other materials provided with the distribution.
+#
+# THIS SOFTWARE IS PROVIDED BY THE AUTHORS ''AS IS'' AND ANY EXPRESS OR IMPLIED
+# WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+# MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
+# NO EVENT SHALL LAB DPR OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
+# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+# POSSIBILITY OF SUCH DAMAGE.
+#
+#
+# The views and conclusions contained in the software and documentation are
+# those of the authors and should not be interpreted as representing
+# official policies, either expressed or implied, of Lab DPR.
+#
+# */
+# /************************************************************************/
+#
+
+arti_path="${LAGO_ARTI}"
+odir=""
+wdir="${PWD}"
+dirlw=0
+energy_bins=20
+distance_bins=20
+altitude=0
+filter=0
+tim=0
+prj=""
+cmd="showers"
+loc=0
+parallel=0
+prims=0
+N=$(/usr/bin/nproc) # number of simultaneous process for paralllel local processing
+N=$((N / 2)) # number of simultaneous process for paralllel local processing
+if [ $N -gt 8 ]; then
+ N=8
+fi
+showhelp() {
+ echo
+ echo -e "$0 version $VERSION"
+ echo
+ echo -e "USAGE $0:"
+ echo
+ echo -e " -o <origin directory> : Origin dir, where the DAT files are located"
+ echo -e " -r <ARTI directory> : ARTI installation directory, generally pointed by \$LAGO_ARTI (default)"
+ echo -e " -w <workding directory> : Working dir, where the analysis will be done (default is current directory, ${wdir})"
+ echo -e " -r <ARTI directory> : ARTI installation directory, generally pointed by \$LAGO_ARTI (default)"
+ echo -e " -e <energy bins> : Number of energy secondary bins (default: $energy_bins"
+ echo -e " -d <distance bins> : Number of distance secondary bins (default: $distance_bins"
+ echo -e " -p <project base name> : Base name for identification of S1 files (don't use spaces). Default: odir basename"
+ echo -e " -k <site altitude, in m> : For curved mode (default), site altitude in m a.s.l. (mandatory)"
+ echo -e " -s <type> : Filter secondaries by type: 1: EM, 2: MU, 3: HD"
+ echo -e " -t <time> : Normalize energy distribution in particles/(m2 s bin), S=1 m2; <t> = flux time (s)."
+ echo -e " -m <bins per decade> : Produce files with the energy distribution of the primary flux per nuclei."
+ echo -e " -j : Produce a batch file for parallel processing. Not compatible with local (-l)"
+ echo -e " -l : Enable parallel execution locally ($N procs). Not compatible with parallel (-j)"
+ echo -e " -? : Shows this help and exit."
+ echo
+}
+echo
+while getopts ':r:w:o:e:p:d:k:s:t:m:lj?' opt; do
+ case $opt in
+ r)
+ arti_path=${OPTARG%/}
+ ;;
+ w)
+ wdir=${OPTARG%/}
+ ;;
+ o)
+ odir=${OPTARG%/}
+ ;;
+ p)
+ prj=$OPTARG
+ ;;
+ e)
+ energy_bins=$OPTARG
+ ;;
+ d)
+ distance_bins=$OPTARG
+ ;;
+ k)
+ altitude=$OPTARG
+ ;;
+ s)
+ filter=$OPTARG
+ ;;
+ m)
+ prims=$OPTARG
+ ;;
+ t)
+ tim=$OPTARG
+ ;;
+ l)
+ loc=1
+ ;;
+ j)
+ parallel=1
+ ;;
+ ?)
+ showhelp
+ exit 1;
+ ;;
+ esac
+done
+
+##################################################
+## YOU SHOULD NOT EDIT ANYTHING BELOW THIS LINE ##
+##################################################
+
+# ERRORS
+
+# is ARTI installed?
+file=$arti_path/analysis/analysis
+if [ ! -f "$file" ]; then
+ echo; echo -e "# ERROR: ARTI analysis executable files not found in $arti_path. Please check and try again"
+ echo
+ showhelp
+ exit 1;
+fi
+
+# are the DAT files in the data directory?
+file=$(ls -1 $odir/DAT??????.bz2 | wc -l)
+if [ ! $file -gt 0 ]; then
+ echo; echo -e "# ERROR: DAT files not found in $odir. Please check and try again"
+ echo
+ showhelp
+ exit 1
+fi
+
+# is working really local? for now just check for /mnt
+if [[ $wdir == *"/mnt/"* ]]; then
+ echo; echo -e "# ERROR: working dir (-w) should be local, not mounted. Please check and try again"
+ echo
+ showhelp
+ exit 1
+fi
+
+# is altitude is given?
+if [ ${altitude} -eq 0 ]; then
+ echo; echo -e "# ERROR: Altitude was not provided. It is mandatory for automatic analysis (curved mode)"
+ echo
+ showhelp
+ exit 1
+fi
+
+# both parallels are not compatible?
+if [ $parallel -gt 0 ] && [ $loc -gt 0 ]; then
+ echo; echo -e "# ERROR: Parallel and local modes are not compatible. Look for -j or -l options."
+ echo
+ showhelp
+ exit 1
+fi
+
+# No errors found, so
+cmd+=" -a $energy_bins"
+cmd+=" -d $distance_bins"
+cmd+=" -c $altitude"
+
+# WARNINGS
+
+if [ $tim -gt 0 ]; then
+ echo; echo -e "# WARNING: Will normalize the secondary flux, S=1 m2; t=$tim"
+ cmd+=" -n 1. $tim"
+fi
+
+if [ $filter -gt 0 ]; then
+ echo; echo -e "# WARNING: Will filter by secondary type $filter (1: EM, 2: MU, 3: HD)"
+ cmd+=" -s $filter"
+fi
+
+if [ $prims -gt 0 ]; then
+ echo; echo -e "# WARNING: Will produce energy histograms of primaries with ${prims} per decade"
+fi
+
+if [ "X$prj" == "X" ]; then
+ prj=$(basename $odir)
+ prj=${prj/S0/S1}
+ echo; echo -e "# WARNING: Project base name not provided. Using $prj"
+fi
+
+cmd+=" $prj"
+
+if [ "X$wdir" == "X$odir" ]; then
+ echo; echo -e "# WARNING: We are running where DAT files are located. Analysis results will be in the same directory."
+ echo -e "# $wdir"
+ dirlw=1
+fi
+
+if [ $loc -gt 0 ]; then
+ if [ $N -gt $file ]; then
+ echo; echo -e "# WARNING: You don't have enough files to analyze in local parallel model (at least $N). Turning off parallel mode."
+ loc=0
+ fi
+fi
+
+## finally...
+wdir=$wdir/$prj
+mkdir $wdir
+cd $wdir
+echo
+echo -e "# Path to ARTI directory = $arti_path"
+echo -e "# Path to DAT files = $odir"
+echo -e "# Path to running directory = $wdir"
+echo -e "# Project base_name = $prj"
+echo -e "# Energy bins = $energy_bins"
+echo -e "# Distance bins = $distance_bins"
+echo -e "# Altitude = $altitude"
+echo -e "# Filtering by type = $filter"
+echo -e "# Normalize flux, S=1 m2; time = $tim"
+echo -e "# Energy bins for primaries = $prims"
+echo -en "# Parallel mode (local) = "
+if [ $loc -gt 0 ]; then
+ echo -e "Local - $N processes"
+else
+ echo -e "Remote"
+fi
+
+# primaries
+for i in ${odir}/DAT??????.bz2; do
+ j=$(basename $i .bz2)
+ u=${j/DAT/}
+ if [ $dirlw -gt 0 ]; then
+ run="bzip2 -d -k $i; echo $j | ${arti_path}/analysis/lagocrkread | ${arti_path}/analysis/analysis -p ${u}; rm ${j}"
+ else
+ run="while ! cp -a $i $wdir/; do sleep 5; done; bzip2 -d $j.bz2; echo $j | ${arti_path}/analysis/lagocrkread | ${arti_path}/analysis/analysis -p ${u}; rm $wdir/${j}"
+ fi
+ echo $run >> $prj.run
+done
+nl=$(cat $prj.run | wc -l)
+if [ $parallel -gt 0 ]; then
+ # parallel mode, just produce the shower analysis file and exit
+ echo "bzcat ${wdir}/*.sec.bz2 | ${arti_path}/analysis/${cmd}" > $prj.shw.run
+ if [ $prims -gt 0 ]; then
+ echo "primaries.sh -w ${wdir} -r ${arti_path} -m ${prims}" > $prj.pri.run
+ fi
+ exit 0
+elif [ $loc -gt 0 ]; then
+ # parallel and local
+ while IFS= read -r line; do
+ ((nr++))
+ ((n=n%N)); ((n++==0)) && wait
+ echo $nr/$nl
+ eval ${line} &>> $nr.log &
+ done < $prj.run
+else
+ # it's local and not parallel
+ while IFS= read -r line; do
+ ((nr++))
+ echo $nr/$nl
+ eval ${line} &>> $nr.log
+ done < $prj.run
+fi
+if [ $loc -gt 0 ]; then
+ echo "Wait for parallel execution termination..."
+ while true; do
+ f=$(find . -iname 'DAT??????' | wc -l)
+ if [ $f -eq 0 ]; then
+ break
+ fi
+ done
+fi
+echo "Shower analysis ..."
+# showers
+echo "Showers: $cmd"
+bzcat ${wdir}/*.sec.bz2 | ${arti_path}/analysis/${cmd}
+
+# primaries histograms
+if [ $prims -gt 0 ]; then
+ echo "Primary analysis ..."
+ primaries.sh -w ${wdir} -r ${arti_path} -m ${prims}
+fi
+
+# final remarks
+# rm $prj.run
+rm *.log
diff --git a/analysis/primaries.sh b/analysis/primaries.sh
new file mode 100755
index 0000000000000000000000000000000000000000..2adf838bf632be43b757d3b7a1c8fccc1203387f
--- /dev/null
+++ b/analysis/primaries.sh
@@ -0,0 +1,141 @@
+#!/bin/bash
+# /************************************************************************/
+# /* */
+# /* Package: ARTI */
+# /* Module: do_showers.sh */
+# /* */
+# /************************************************************************/
+# /* Authors: Hernán Asorey */
+# /* e-mail: hernanasorey@cnea.gob.ar */
+# /* */
+# /************************************************************************/
+# /* Comments: Script to automatize the analysis of simulated showers */
+# /* */
+# /************************************************************************/
+# /*
+#
+#
+# Copyright 2021
+# Hernán Asorey
+# ITeDA (CNEA), Argentina
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are
+# met:
+#
+# 1. Redistributions of source code must retain the above copyright
+# notice, this list of conditions and the following disclaimer.
+#
+# 2. Redistributions in binary form must reproduce the above copyright
+# notice, this list of conditions and the following disclaimer in the
+# documentation and/or other materials provided with the distribution.
+#
+# THIS SOFTWARE IS PROVIDED BY THE AUTHORS ''AS IS'' AND ANY EXPRESS OR IMPLIED
+# WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+# MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
+# NO EVENT SHALL LAB DPR OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
+# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+# HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+# STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+# POSSIBILITY OF SUCH DAMAGE.
+#
+#
+# The views and conclusions contained in the software and documentation are
+# those of the authors and should not be interpreted as representing
+# official policies, either expressed or implied, of Lab DPR.
+#
+# */
+# /************************************************************************/
+#
+
+VERSION="v1r0" # First release, mar 20 abr 2021 11:00:00 CEST
+
+wdir=${PWD}
+arti_path=${LAGO_ARTI}
+prj=""
+prims=10
+
+showhelp() {
+ echo
+ echo -e "$0 version $VERSION"
+ echo
+ echo -e "USAGE $0:"
+ echo
+ echo -e " -w <working directory> : Working directory, where the pri.bz2 iles are located"
+ echo -e " -r <ARTI directory> : ARTI installation directory, generally pointed by \$LAGO_ARTI (default)"
+ echo -e " -m <bins per decade> : Produce files with the energy distribution of the primary flux per nuclei. Not compatible with parallel"
+ echo -e " -? : Shows this help and exit."
+ echo
+}
+echo
+while getopts 'w:r:p:m:?' opt; do
+ case $opt in
+ w)
+ wdir=$OPTARG
+ ;;
+ r)
+ arti_path=$OPTARG
+ ;;
+ m)
+ prims=$OPTARG
+ ;;
+ ?)
+ showhelp
+ exit 1;
+ ;;
+ esac
+done
+
+##################################################
+## YOU SHOULD NOT EDIT ANYTHING BELOW THIS LINE ##
+##################################################
+
+# ERRORS
+file=$(ls -1 $wdir/*.pri.bz2 | head -1)
+echo $file
+if [ ! -f "$file" ]; then
+ echo; echo -e "# ERROR: pri.bz2 files not found in $wdir. Please check and try again"
+ echo
+ showhelp
+ exit 1
+fi
+
+file=$arti_path/analysis/analysis
+if [ ! -f "$file" ]; then
+ echo; echo -e "# ERROR: ARTI analysis executable files not found in $arti_path. Please check and try again"
+ echo
+ showhelp
+ exit 1;
+fi
+
+# WARNINGS
+
+pass=1
+if [ "X$PWD" == "X$wdir" ]; then
+ pass=0
+else
+ echo; echo -e "# WARNING: Not running where pri.bz2 files are located. At the end will move all files to $wdir"
+fi
+
+## finally...
+echo
+echo -e "# Path to ARTI directory = $arti_path"
+echo -e "# Path to DAT files = $wdir"
+echo -e "# Energy bins for primaries = $prims"
+
+# primaries histograms
+ids=$(for i in *.pri.bz2; do j=${i/.pri.bz2/}; k=${j:2}; echo $k; done | sort | uniq)
+echo ${ids}
+echo
+for i in $ids; do
+ echo -n "${i} "
+ bzcat ??${i}.pri.bz2 | grep -v "#" | awk '{print log($2)/log(10.)}' | sort -g | awk -v bins=${prims} -v id=${i} 'BEGIN{n=0; mine=100000; maxe=-100000; bins = bins * 1.}{t[int($1*bins)]++; n++; if ($1 < mine) mine=$1; if ($1 > maxe) maxe=$1;}END{printf("# # # prt\n");printf("# # Primary energy histogram for %06d using %d bins per decade\n", id, bins);printf("# # Three column format is:\n# # energy_bin total_per_bin fraction_per_bin\n"); for (i in t) {print 10**(i/bins), t[i], t[i]*1./n; frc+=t[i]*1./n;} printf("# # Total primaries: %ld (%.2f) Emin=%.2f GeV; Emax=%.2f GeV\n", n, frc, 10**mine, 10**maxe);}' > "00${i}.prt"
+done
+echo
+if [ $pass -gt 0 ]; then
+ mv *.prt $wdir
+fi
diff --git a/analysis/showers.cc b/analysis/showers.cc
index 53b41c10ea64b5241cb5e286aa4631bb5f9e7286..58ca34c29250460727df3c760a7d82fdaf4bb469 100644
--- a/analysis/showers.cc
+++ b/analysis/showers.cc
@@ -65,13 +65,14 @@ using namespace std;
int iverbose = 0, iforce = 0, ifilter = 0, icurve = 0, idistance=0, ianalysis=0, inorm = 0, igeo = 0;
int particle=0;
bool print=true;
-FILE *fi, *shw, *hst, *dst;
+FILE *fi, *shw, *hst, *dst, *dse;
double x, y, z, h, hInM;
double r_earth=637131500.;
double GeV2keV=1.0e6;
+double kev2GeV=1.0e-6;
double maxPInGeV = 1.0e7; //10^16 eV is the upper limit for p particles (not expected)
double cm2m=0.01;
-double maxDInM = 6e4; // 60 km should be enough
+double maxDInM = 5e4; // 50 km should be enough
int cGeoCdef=3;
int cGeoC = cGeoCdef; // column from geomagnetic cutoff file
@@ -84,6 +85,49 @@ int masaI[masaN] = {1, 14, 402, 703, 904, 1105, 1206, 1407, 1608, 1909, 2010, 23
double masas[masaN] = {0., 0.938272, 3.73338, 6.53478, 8.39429, 10.25361, 11.17879, 13.04507, 14.89833, 17.68991, 18.61736, 21.40802, 22.33628, 25.12634, 26.05532, 28.84497, 29.77460, 32.56398, 37.21075, 36.28343, 37.21424, 41.86053, 44.64837, 47.44020, 48.36813, 51.15981, 52.08852};
+double mass(int id) {
+// return mass in GeV for the Corsika id secondary
+ switch (id) {
+ //1=foton, 2=e+, 3=e-, 5=mu+, 6=mu-, 7=pi0, 8=pi+, 9=pi-, 13=n, 14=p, 15=bar-p, other
+ case 1:
+ return (0.); //photon
+ break;
+ case 2:
+ case 3:
+ return (511.*kev2GeV);
+ break;
+ case 5:
+ case 6:
+ return (0.10566);
+ break;
+ case 7:
+ return (0.13498);
+ break;
+ case 8:
+ case 9:
+ return (0.13957);
+ break;
+ case 13:
+ case 25:
+ return (0.93957);
+ break;
+ case 14:
+ case 15:
+ return (0.938272);
+ break;
+ default:
+ if (id < 100)
+ return 1.;
+ else {
+ int idm = 0;
+ for (int i=0; i<masaN; i++)
+ if (masaI[i] == id)
+ idm = i;
+ return masas[idm];
+ }
+ return 1; // just in case
+ }
+}
int Open(char *nfi) {
char tmpc[256];
@@ -111,8 +155,11 @@ void Usage(char *prog, int iverbose=0) {
cout << " -s <type> : filter secondaries by type: 1: EM, 2: MU, 3: HD"<<endl;
cout << " -a <n> : Energy distribution of secondaries per type, with <n> bins per decade."<<endl;
cout << " -d <n> : Distance distribution of secondaries per type, with <n> bins per decade."<<endl;
+ cout << " : If normalization is also selected (-n) the distribution is weighted by the ring area."<<endl;
+ cout << " : In that case, produce also a new file (dse) with the energy density distribution at ground."<<endl;
cout << " -n <a> <t> : Normalize energy distribution in particles/(m2 s bin), <a>=detector area (m2), <t> = flux time (s)." <<endl;
- cout << " -c <obs_lev in m>: Enable curved mode: x'y' will be converted to local coordinates xy"<<endl;
+ cout << " : Change distance distribution by normalizing by the increasing area of the ring (r, r+dr)." <<endl;
+ cout << " -c <obs_lev (m)>: Enable curved mode: x'y' will be converted to local coordinates xy"<<endl;
cout << " Observation level (<obs_lev>) should be given in m a.s.l." << endl;
cout << " -g <file> <col> : Include geomagentic effects. Read rigidities from column <R> of <file>. Default R= " << cGeoCdef <<endl;
cout << " 3=R_U, 4=R_C 5=R_L"<<endl;
@@ -135,10 +182,16 @@ void curved(double xp, double yp) {
y = d*sin(f)*cm2m;
}
-inline double momemtum(double px, double py, double pz) {
+inline double momentum(double px, double py, double pz) {
return (GeV2keV*sqrt(px*px+py*py+pz*pz));
}
+inline double energy(double m, double pkeV) {
+// Energy in GeV for a particle of mass m and momentum m
+ double p = pkeV * kev2GeV; // Asuming pkeV comes from momentum function (in keV)
+ return (sqrt(m * m + p * p));
+}
+
inline double distance(double r1, double r2, double r3) {
return (sqrt(r1*r1+r2*r2+r3*r3));
}
@@ -147,6 +200,14 @@ inline double log10(double x) {
return (log(x)/log(10.));
}
+inline double ring(double r1, double r2) {
+ if (r2 > r1)
+ return M_PI*(r2*r2 - r1*r1);
+ else
+ return M_PI*(r1*r1 - r2*r2);
+}
+
+
int main (int argc, char *argv[]) {
char nfi[256]; // filename container
@@ -312,8 +373,14 @@ int main (int argc, char *argv[]) {
fprintf(stderr,"Failed to open dst (distance histogram) file. Abort.\n");
exit(1);
}
+ snprintf(nfi, 256, "%s.dse", ifile);
+ if (inorm) {
+ if ((dse = fopen(nfi,"w"))==NULL) {
+ fprintf(stderr,"Failed to open dse (distance energy histogram) file. Abort.\n");
+ exit(1);
+ }
+ }
}
-
fprintf(shw, "# # # shw\n");
if (icurve)
fprintf(shw, "# # CURVED mode is ENABLED and observation level is %.0lf m a.s.l.\n", hInM);
@@ -345,6 +412,7 @@ int main (int argc, char *argv[]) {
//1=foton, 2=e+, 3=e-, 5=mu+, 6=mu-, 7=pi0, 8=pi+, 9=pi-, 13=n, 14=p, 15=bar-p, other
}
if (idistance) {
+ // DST
fprintf(dst, "# # # dst\n");
if (iforce)
fprintf(dst, "# # # WARNING: force mode is enable. Analysis done in FLAT mode.\n");
@@ -353,12 +421,30 @@ int main (int argc, char *argv[]) {
if (igeo)
fprintf(dst, "# # Geomagnetic effects were included: %d values were read from file %s\n", iGeoN, nfg);
fprintf(dst, "# # This is the Histogram of secondary distance file - ARTI %s\n", VERSION);
- fprintf(dst, "# # Logaritmic distance scale. Resolution used: %d bins per distance decade\n", int(resolution));
+ fprintf(dst, "# # Logarithmic distance scale. Resolution used: %d bins per distance decade\n", int(resolution));
if (inorm)
- fprintf(dst, "# # Number of particles are divided by detector area (%.4f m^2) and flux time (%.2f s)\n", area, fluxTime);
+ fprintf(dst, "# # Number of particles are divided by detector area (%.4f m^2), ring (r,r+dr) area and flux time (%.2f s).\n", area, fluxTime);
fprintf(dst, "# # 14 column format is:\n");
fprintf(dst, "# # distance_in_bin(m) N_phot N_e+ N_e- N_mu+ N_mu- N_pi0 N_pi+ N_pi- N_n N_p N_pbar N_others Total_per_bin\n");
//1=foton, 2=e+, 3=e-, 5=mu+, 6=mu-, 7=pi0, 8=pi+, 9=pi-, 13=n, 14=p, 15=bar-p, other
+
+ // DSE
+ if (inorm) {
+ fprintf(dse, "# # # dse\n");
+ if (iforce)
+ fprintf(dse, "# # # WARNING: force mode is enable. Analysis done in FLAT mode.\n");
+ if (icurve)
+ fprintf(dse, "# # CURVED mode is ENABLED and observation level is %.0lf m a.s.l.\n", hInM);
+ if (igeo)
+ fprintf(dse, "# # Geomagnetic effects were included: %d values were read from file %s\n", iGeoN, nfg);
+ fprintf(dse, "# # This is the Histogram of secondary energy distance distribution file - ARTI %s\n", VERSION);
+ fprintf(dse, "# # Logarithmic distance scale. Resolution used: %d bins per distance decade\n", int(resolution));
+ if (inorm)
+ fprintf(dse, "# # Total energy (in GeV) of particles at each bin are divided by detector area (%.4f m^2), ring area (r,r+dr) and flux time (%.2f s)\n", area, fluxTime);
+ fprintf(dse, "# # 14 column format is:\n");
+ fprintf(dse, "# # distance_in_bin(m) E_phot E_e+ E_e- E_mu+ E_mu- E_pi0 E_pi+ E_pi- E_n E_p E_pbar E_others Total_E_per_bin\n");
+ //1=foton, 2=e+, 3=e-, 5=mu+, 6=mu-, 7=pi0, 8=pi+, 9=pi-, 13=n, 14=p, 15=bar-p, other
+ }
}
if (iverbose) {
fprintf (stderr, "Verbosity level set to: %d.\n", iverbose);
@@ -380,10 +466,13 @@ int main (int argc, char *argv[]) {
nlogd = int(log10(maxDInM)*resdist);
long int histod[nlogd][nbin];
long int nld[nlogd];
+ double histode[nlogd][nbin];
+ double nlde[nlogd];
long int maxDerr = 0, maxEerr = 0;
int prmId, prmX, prmZ, prmA;
double prmEn, prmTh, prmPh;
double prmMa, prmIm;
+ double secp = 0.;
bool prmGeo = true;
long int geoDisShw = 0, geoDisSec = 0;
long int totpart=0, totbin = 0;
@@ -396,13 +485,19 @@ int main (int argc, char *argv[]) {
for (int j=0; j<nlogd; j++) {
histod[j][i] = 0;
}
+ for (int j=0; j<nlogd; j++) {
+ histode[j][i] = 0.;
+ }
+
nid[i] = 0;
}
for (int j=0; j<nloge; j++)
nle[j] = 0;
for (int j=0; j<nlogd; j++)
nld[j] = 0;
-
+ for (int j=0; j<nlogd; j++)
+ nlde[j] = 0;
+
long int lines = 0, ps = 0, shower_id = 0;
while(fgets(line,250,fi)) {
if (!(++lines % 10000))
@@ -450,6 +545,7 @@ int main (int argc, char *argv[]) {
if (prmId==masaI[i])
prmX = i;
if (prmX < 0) {
+ fprintf(stderr, " Error: Can't find primary. PrmId: %04d. Please check. \nLine: %s\n", prmId, line);
fprintf(stderr, " Error: Can't find primary. PrmId: %04d. Please check. \nLine: %s\n", prmId, line);
exit(1);
}
@@ -547,14 +643,15 @@ int main (int argc, char *argv[]) {
ibin=11; //other hadrons
break;
}
+ secp = momentum(d[1], d[2], d[3]);
if (ianalysis) {
- lp = log10(momemtum(d[1], d[2], d[3]));
+ lp = log10(secp);
iloge=int(lp*resolution);
if (iloge > nloge) {
iloge=nloge;
maxEerr++;
if (iverbose>1)
- cerr << "Warning! Particle momemtum is " << momemtum(d[1], d[2], d[3]) << "which is bigger than " << maxPInGeV << " GeV. Should check. \nLine: " << line << endl;
+ cerr << "Warning! Particle momentum is " << secp << "which is bigger than " << maxPInGeV << " GeV. Should check. \nLine: " << line << endl;
}
if (iloge < minbine)
minbine=iloge;
@@ -567,8 +664,6 @@ int main (int argc, char *argv[]) {
ld = log10(distance(x, y, (z-hInM)));
if (ld < 0.)
ld = 0.;
-// cout
-// << distance(x, y, (z-hInM)) << endl;
ilogd=int(ld*resdist);
if (ilogd > nlogd) {
ilogd=nlogd;
@@ -582,6 +677,8 @@ int main (int argc, char *argv[]) {
maxbind=ilogd;
histod[ilogd][ibin]++;
nld[ilogd]++;
+ histode[ilogd][ibin] += energy(mass(id), secp);
+ nlde[ilogd] += energy(mass(id), secp);
}
nid[ibin]++;
totbin++;
@@ -597,9 +694,10 @@ int main (int argc, char *argv[]) {
if (igeo)
fprintf(stderr,"Including geomagnetic effects, %ld secondaries from %ld showers were discarded.\n", geoDisSec, geoDisShw);
- double p=0.,r=0;
+ double p = 0., r = 0., r1 = 0., rarea = 0.;
+
double norm=area*fluxTime;
-
+ // energy
if (ianalysis) {
for (int i=minbine; i<=maxbine; i++) {
p=pow(10., (i/resolution))/GeV2keV;
@@ -637,18 +735,27 @@ int main (int argc, char *argv[]) {
fprintf(hst, "# Totals: EM:MU:NE:HD= %ld:%ld:%ld:%ld\n", (nid[0]+nid[1]+nid[2]), (nid[3]+nid[4]), (nid[8]), (nid[5]+nid[6]+nid[7]+nid[9]+nid[10]+nid[11]));
fprintf(hst, "# Ratios: EM:MU:NE:HD= %02.3f:%02.3f:%02.3f:%02.3f\n", (nid[0]+nid[1]+nid[2])/tot, (nid[3]+nid[4])/tot, (nid[8])/tot, (nid[5]+nid[6]+nid[7]+nid[9]+nid[10]+nid[11])/tot);
}
+ // warning, normalization depends also on distance due to the ring area...
if (idistance) {
for (int i=minbind; i<=maxbind; i++) {
r=pow(10., (i/resdist));
+ r1=pow(10., ((i+1)/resdist));
+ rarea = norm * ring(r1, r);
fprintf(dst, "%.6e ", r);
+ if (inorm)
+ fprintf(dse, "%.6e ", r);
for (int j=0; j<nbin; j++) {
- if (inorm)
- fprintf(dst, "%.7e ", histod[i][j]/norm);
+ if (inorm) {
+ fprintf(dst, "%.7e ", histod[i][j]/rarea);
+ fprintf(dse, "%.7e ", histode[i][j]/rarea);
+ }
else
fprintf(dst, "%ld ", histod[i][j]);
}
- if (inorm)
- fprintf(dst, "%.7e\n", nld[i]/norm);
+ if (inorm) {
+ fprintf(dst, "%.7e\n", nld[i]/rarea);
+ fprintf(dse, "%.7e\n", nlde[i]/rarea);
+ }
else
fprintf(dst, "%ld\n", nld[i]);
}
@@ -659,7 +766,7 @@ int main (int argc, char *argv[]) {
if (inorm) {
fprintf(dst, "#TOT_NORM: ");
for (int j=0; j<nbin; j++)
- fprintf(dst, "%.7e ", nid[j]/norm);
+ fprintf(dst, "%.7e ", nid[j]/rarea);
fprintf(dst, "\n");
}
double tot = totbin * 1.0;
@@ -669,7 +776,7 @@ int main (int argc, char *argv[]) {
fprintf(dst, "\n");
fprintf(dst, "# Total number of binned particles: %ld", totbin);
if (inorm)
- fprintf(dst, " (corresponding to $%.2f$\\,particles\\, m$^{-2}$\\,s$^{-1}$.)", tot/norm);
+ fprintf(dst, " (corresponding to $%.2f$\\,particles\\,s$^{-1}$, distributed in a total area of %.0f\\,m$^{-2}$)", tot/norm, ring(r, 0.));
fprintf(dst, "\n");
fprintf(dst, "# Totals: EM:MU:NE:HD= %ld:%ld:%ld:%ld\n", (nid[0]+nid[1]+nid[2]), (nid[3]+nid[4]), (nid[8]), (nid[5]+nid[6]+nid[7]+nid[9]+nid[10]+nid[11]));
fprintf(dst, "# Ratios: EM:MU:NE:HD= %02.3f:%02.3f:%02.3f:%02.3f\n", (nid[0]+nid[1]+nid[2])/tot, (nid[3]+nid[4])/tot, (nid[8])/tot, (nid[5]+nid[6]+nid[7]+nid[9]+nid[10]+nid[11])/tot);
@@ -686,8 +793,9 @@ int main (int argc, char *argv[]) {
}
if (idistance) {
if (maxDerr)
- fprintf(dst, "# # WARNING: %ld particles (%.1f%%) overpass max distance bin (located at last bin)\n", maxDerr, (100.*maxDerr)/totpart);
+ fprintf(dst, "# # WARNING: %ld particles (%.1f%%) overpass max distance bin (%d km), and were added to last bin.\n", maxDerr, (100.*maxDerr)/totpart, int(maxDInM/1000.));
fprintf(dst, "# # Selected particles: %ld (binned, %ld), showers: %ld, lines: %ld\n", totpart, totbin, shower_id, lines);
fclose(dst);
+ fclose(dse);
}
}
diff --git a/sims/generate_spectra.pl b/sims/generate_spectra.pl
index 76a6ebbf69d866d93870257622f2a3221de3b6f1..0fde837460ce0108226160b39a335e7f50a37f3d 100755
--- a/sims/generate_spectra.pl
+++ b/sims/generate_spectra.pl
@@ -51,7 +51,7 @@
# */
# /************************************************************************/
#
-$VERSION="v1r1";
+$VERSION="v1r2";
# defaults
use Switch;
@@ -140,7 +140,7 @@ Optional:
-u <user name> For CORSIKA simulation (Default: none)
-t <time in sec> Flux time ( Default: 3600s)
-s <site> Choice predefined site for simulation (Default: unknow)
- - Predefined sites: hess|sac|etn|ber|bga|lim|glr|mch|mge|and|mpc|cha|cid|mor|ccs|lsc|mbo)
+ - Predefined sites: check the code)
- Predefined parameters: altitude, BX, BZ, and Atmospheric Model.
-k <altitude> Fix altitude even for predefined sites (It cannot be 0)
-c <modatm> Fix Atmospheric Model even for predefined sites.
@@ -291,124 +291,239 @@ $userllimit=$llimit;
#all predefined sites seems to be ATMOSPHERE as default mode
$atmcrd = "ATMOSPHERE";
+
+# New set of LAGO sites in preparation for EOSC challenge.
+# For pre-challenge we will simulate some of this new sites.
+# HA - Apr 08, 2021
+
switch ($site) {
- case "hess" {
- $modatm="E10";
- $altitude=1800e2;
- $bx=12.5;
- $bz=-25.9;
+ case "sawb" {
+ $modatm="E5";
+ $altitude=20000;
+ $bx=19.366;
+ $bz=-30.222;
+ }
+ case "mapi" {
+ $modatm="E5";
+ $altitude=1000;
+ $bx=19.309;
+ $bz=-28.693;
}
- case "sac" {
- #$modatm="E32"; ### OJO ajrm no funciona "ATMOSPHERE E32" , de hecho ATMOSPHERE no parece que le guste mucho
- $modatm="10"; ### OJO ajrm si pongo esto funciona y me lo cambia a "ATMOD 10" no sé donde
- $altitude=3700e2;
- $bx=20.94;
- $bz=-8.91;
- }
- case "etn" {
- $modatm="E2";
- $altitude=3000e2;
- $bx=27.7623;
- $bz=36.0667;
+ case "brc" {
+ $modatm="E3";
+ $altitude=86500;
+ $bx=18.952;
+ $bz=-17.05;
+ }
+ case "bue" {
+ $modatm="E2";
+ $altitude=1000;
+ $bx=17.09;
+ $bz=-14.673;
+ }
+ case "kna" {
+ $modatm="E2";
+ $altitude=34700;
+ $bx=19.56;
+ $bz=-13.333;
}
- case "ber" {
- $modatm="E1";
- $altitude=3450e2;
- $bx=26.9814;
- $bz=17.1054;
+ case "lsc" {
+ $modatm="E2";
+ $altitude=2800;
+ $bx=19.975;
+ $bz=-11.826;
+ }
+ case "tuc" {
+ $modatm="E2";
+ $altitude=43000;
+ $bx=19.487;
+ $bz=-10.647;
+ }
+ case "asu" {
+ $modatm="E1";
+ $altitude=13600;
+ $bx=18.233;
+ $bz=-11.73;
+ }
+ case "sao" {
+ $modatm="E1";
+ $altitude=76000;
+ $bx=16.484;
+ $bz=-14.48;
+ }
+ case "vcp" {
+ $modatm="E1";
+ $altitude=64000;
+ $bx=16.751;
+ $bz=-14.104;
+ }
+ case "lpb" {
+ $modatm="E2";
+ $altitude=363000;
+ $bx=22.362;
+ $bz=-4.708;
}
- case "lim" {
- $modatm="E2";
- $altitude=168e2;
- $bx=25.28;
- $bz=-0.046;
- }
- case "glr" {
- $modatm="E1";
- $altitude=4276e2;
- $bx=27.0750;
- $bz=11.7728;
- $arrang="0";
- }
- case "mch" {
- $modatm="E1";
- $altitude=2650e2;
- $bx=27.1762;
- $bz=14.6184;
- $arrang="0";
+ case "cha" {
+ $modatm="E2";
+ $altitude=523000;
+ $bx=22.38;
+ $bz=-4.62;
}
- case "bga" {
- $modatm="E1";
- $altitude=950e2;
- $bx=27.0263;
- $bz=17.1760;
+ case "cuz" {
+ $modatm="E1";
+ $altitude=340000;
+ $bx=23.676;
+ $bz=-2.02;
}
- case "mge" {
- $modatm="19";
- $altitude=1400e2;
- $bx=20.4367;
- $bz=-11.8217;
+ case "lim" {
+ $modatm="E1";
+ $altitude=16800;
+ $bx=24.723;
+ $bz=-0.476;
+ }
+ case "cpv" {
+ $modatm="E1";
+ $altitude=55000;
+ $bx=21.174;
+ $bz=-12.612;
+ }
+ case "serb" {
+ $modatm="E1";
+ $altitude=275000;
+ $bx=26.577;
+ $bz=8.499;
+ }
+ case "QUIE" {
+ $modatm="E1";
+ $altitude=285000;
+ $bx=26.717;
+ $bz=9.971;
}
- case "brc" {
- $modatm="E3";
- $altitude=800e2;
- $bx=19.234;
- $bz=-17.068;
+ case "bga" {
+ $modatm="E1";
+ $altitude=95000;
+ $bx=26.793;
+ $bz=16.055;
+ }
+ case "pam" {
+ $modatm="E1";
+ $altitude=234200;
+ $bx=26.743;
+ $bz=15.955;
+ }
+ case "gua" {
+ $modatm="E1";
+ $altitude=149000;
+ $bx=27.532;
+ $bz=24.932;
+ }
+ case "tac" {
+ $modatm="E1";
+ $altitude=406000;
+ $bx=27.53;
+ $bz=25.353;
+ }
+ case "chia" {
+ $modatm="E2";
+ $altitude=52200;
+ $bx=27.405;
+ $bz=27.024;
+ }
+ case "sng" {
+ $modatm="E2";
+ $altitude=455000;
+ $bx=27.358;
+ $bz=28.038;
}
case "and" {
- $modatm="19";
- $altitude=4200e2;
- $bx=19.6922;
- $bz=-14.2420;
- }
- case "mpc" {
- # Marcapomacocha
- $modatm="E1";
- $altitude=4500e2;
- $bx=24.9599;
- $bz=0.4124;
+ $modatm="E2";
+ $altitude=420000;
+ $bx=19.658;
+ $bz=-11.951;
}
- case "cha" {
- # Chacaltaya
- $modatm="E2";
- $altitude=5230e2;
- $bx=23.0386;
- $bz=-3.9734;
- }
- case "cid" {
- # CIDA
- $modatm="E1";
- $altitude=3600e2;
- $bx=26.8464;
- $bz=+18.1604;
- }
- case "mor" {
- # Mordor
- $modatm="E1";
- $altitude=4400e2;
- $bx=26.8340;
- $bz=+18.2004;
+ case "mge" {
+ $modatm="E2";
+ $altitude=145000;
+ $bx=18.987;
+ $bz=-14.326;
}
- case "lsc" {
- # La Serena
- $modatm="E2";
- $altitude=28e2;
- $bx=20.29;
- $bz=-11.74;
- }
- case "mbo" {
- # Base Marambio
- $modatm="E5";
- $altitude=196e2;
- $bx=19.6571;
- $bz=-30.5809;
- }
- case "ccs" {
- #Caracas, data provided by Jose Antonio López, UCV, 10.486004N -66.894461W
- $modatm="E1";
- $altitude=900E2;
- $bx=26.7364;
- $bz=+18.6777;
+ case "ber" {
+ $modatm="E1";
+ $altitude=345000;
+ $bx=26.751;
+ $bz=16.029;
+ }
+ case "sac" {
+ $modatm="E2";
+ $altitude=482000;
+ $bx=20.06;
+ $bz=-9.616;
+ }
+ case "cop" {
+ $modatm="E2";
+ $altitude=200000;
+ $bx=19.04;
+ $bz=-15.493;
+ }
+ case "sgd" {
+ $modatm="E2";
+ $altitude=25000;
+ $bx=18.108;
+ $bz=-16.891;
+ }
+ case "casp" {
+ $modatm="E2";
+ $altitude=260000;
+ $bx=19.473;
+ $bz=-12.436;
+ }
+ case "ppet" {
+ $modatm="E2";
+ $altitude=350000;
+ $bx=19.139;
+ $bz=-14.27;
+ }
+ case "mad" {
+ $modatm="E2";
+ $altitude=70000;
+ $bx=25.647;
+ $bz=36.933;
+ }
+ case "truj" {
+ $modatm="E2";
+ $altitude=56000;
+ $bx=26.223;
+ $bz=35.844;
+ }
+ case "pozn" {
+ $modatm="E2";
+ $altitude=10000;
+ $bx=18.598;
+ $bz=46.45;
+ }
+ case "juli" {
+ $modatm="E2";
+ $altitude=10000;
+ $bx=19.676;
+ $bz=44.928;
+ }
+ case "sudb" {
+ $modatm="E2";
+ $altitude=210000;
+ $bx=17.037;
+ $bz=51.991;
+ }
+ case "viri" {
+ $modatm="E2";
+ $altitude=285000;
+ $bx=19.061;
+ $bz=-16.289;
+ }
+ case "ima" {
+ $modatm="E1";
+ $altitude=460000;
+ $bx=22.969;
+ $bz=-3.79;
}
}#switch
@@ -419,18 +534,18 @@ unless ($ifixmodatm) {
unless ($ifixalt) {
while (!$altitude) {
$altitude = get("Observation level above sea level [cm]",0,"OBSLEV", $altitude);
- if (!$altitude) { print STDERR "ERROR: Observation level is mandatory > 0\n"; }
+ if (!$altitude) { die "ERROR: Observation level is mandatory > 0\n"; }
}
}
while (!$bx) {
$bx=get("Horizontal comp. of the Earth's mag. field (MAGNET) [North,muT],\nsee values at http://www.ngdc.noaa.gov/geomagmodels/struts/calcIGRFWMM",0,"BX", $bx);
- if (!$bx) { print STDERR "ERROR: BX is mandatory\n"; }
+ if (!$bx) { die "ERROR: BX is mandatory\n"; }
}
while (!$bz) {
$bz=get("Vertical comp. of the Earth's mag. field (MAGNET) [downwards,muT]",0,"BZ", $bz);
- if (!$bz) { print STDERR "ERROR: BZ is mandatory\n"; }
+ if (!$bz) { die "ERROR: BZ is mandatory\n"; }
}
if ($ifixalt && $fixalt) {
@@ -439,7 +554,7 @@ if ($ifixalt && $fixalt) {
if ($ifixmodatm && $fixmodatm) {
$modatm=$fixmodatm;
}
-$modatm=uc($modatm);
+# $modatm=uc($modatm);
# just in case :)
######################################################
@@ -534,7 +649,7 @@ $bx
$bz
F
F
-T
+F
F
F", $id, $llimit, $ulimit);
close($fi);
@@ -562,4 +677,4 @@ foreach $z (sort {$b <=> $a} keys %spc) {
}
close($fh);
print"\n";
-system ("clear; echo 'Fluxes'; cat $file");
+system ("echo; echo; echo 'Fluxes'; cat $file");
diff --git a/sims/rain.pl b/sims/rain.pl
index b7cf287e125afb0dbb50a7f4afbd8f436df0bd9a..c15bd172ac414381d47cef0809da12cd28044a28 100755
--- a/sims/rain.pl
+++ b/sims/rain.pl
@@ -60,7 +60,7 @@ $tmp = "";
$batch = 0;
$runmode = 0;
$wdir = "x";
-$crk_ver = "75600";
+$crk_ver = "77402";
$heim = "QGSII";
$debug = 0;
$help = 0;
@@ -81,6 +81,7 @@ $cherenkov = 0;
$grid = 0;
$imuaddi = 0;
$nofruns = 1;
+$ecutshe = 800.;
sub get {
my $question = $_[0];
@@ -192,6 +193,7 @@ while ($_ = $ARGV[0]) {
}
$package="corsika".$crk_ver."Linux_".$heim."_gheisha";
+
if ($ithin) {
$package= $package . "_thin";
}
@@ -574,20 +576,25 @@ for ($i=0; $i<$nofruns; $i++) {
# using external atmospheres bernlhor
$atmcrd = "ATMOSPHERE";
$modatm =~ s/E//g;
- $modatm .= " Y";
+ $modatm .= " Y";
+ } else {
+ if (uc(substr($modatm,0,4)) eq "GDAS") {
+ # gdas model
+ $atmcrd = "ATMFILE";
+ $modatm = "'" . lc($modatm) . "'";
+ }
}
#LAGO ECUTS
# @ecuts=(0.05,0.05,1E-4,1E-4);
- @ecuts=(0.05, 0.05, 0.00005, 0.00005);
+# @ecuts=(0.05, 0.05, 0.00005, 0.00005);
+ @ecuts=(0.05, 0.01, 0.00005, 0.00005);
if ($highsec) {
- @ecuts=(10., 10., 10., 10.); # using 10 GeV
- # if you want to use your own cuts, please add a site and use an if like
- # the one show below for andes:
- if ($site eq "and") {
- @ecuts=(800.,800.,800.,800.); # ecuts for ANDES
- }
- }
+ @ecuts=($ecutshe, $ecutshe, $ecutshe, $ecutshe);
+ if ($elow < $ecutshe) {
+ $elow = $ecutshe;
+ }
+}
#true of false
unless ($batch) {
@@ -623,6 +630,8 @@ for ($i=0; $i<$nofruns; $i++) {
$cards="";
$plotshs="";
$direct2 = $direct;
+ $llongi = "F";
+ $datbas = "F";
if ($grid) {
$direct2 = "."
}
@@ -670,7 +679,7 @@ ECUTS $ecuts[0] $ecuts[1] $ecuts[2] $ecuts[3]
$muadditxt
MUMULT T
-MAXPRT 0
+MAXPRT 1
ELMFLG F T
LONGI $llongi 20. T T
ECTMAP 1.E3
@@ -688,7 +697,7 @@ EXIT
";
}
else {
- $cards="RUNNR $runnr
+ $cards="RUNNR $runnr
EVTNR $evtnr
NSHOW $nshow
@@ -711,10 +720,10 @@ ECUTS $ecuts[0] $ecuts[1] $ecuts[2] $ecuts[3]
$curvout
$muadditxt
MUMULT T
-MAXPRT 0
+MAXPRT 1
ELMFLG F T
LONGI $llongi 10. T T
-ECTMAP 1.E3
+ECTMAP 1.E11
$plotshs
DIRECT $direct2/