4. Subroutines and variables

4.1. User accessible types

type  coord

During the paticle tracking, this system is used.

Type fields:
  • % r (3) [real*8]

  • % sys [character(4)] :: which system. one of ‘xyz’, ‘llh’ or, ‘sph’.

Header :

Zcoord.h

type  position

location of a particle

Type fields:
  • % xyz [coord] :: in xyz

  • % radiallen [real(8)] :: in m . radial length from the center of Earth

  • % depth [real(8)] :: in kg/m2 depth

  • % height [real(8)] :: in m. vertical height(from sea level

  • % colheight [real(8)] :: in m. where the latest nuclear collision took place. (iniitial value is very large value).

Header :

Zpos.h

type  direc
Type fields:
  • % w [coord]

  • % coszenith [real(8)] :: cos of the zenith angle. it is defined as follows: Let’s assume w and position are given in xyz sytem.

Header :

Zdirec.h

type  fmom
Type fields:
  • % p (4) [real(8)] :: four momentum in GeV. p(1) is x component. Note. Momentum is given in the Earth xyz system.

type  ptcl

particle at production

Type fields:
  • % fm [fmom] :: 4 momentum

  • % mass [real(8)] :: mass

  • % code [integer(2)] :: particle code

  • % subcode [integer(2)] :: used mainly to identify paticle/antiparticle if the difference is important. To set particle, “ptcl” is used. anti-partilce, ‘antip” is used for particles For particles of which partilce/antiparticle nature can be judded by its code and charge, the user need not specify it when using cmkptc subroutine. give 0. subcode for gamma ray may be used to identify brems gamma and direct gamma by kdiretg, kcasg

  • % charge [integer(2)] :: charge

Header :

Zptcl.h

type  track

full particle attributes in Cosmos

Type fields:
  • % p [ptcl] :: basic ptcl attributes.

  • % pos [position] :: position

  • % t [real(8)] :: time in length/beta (m)

  • % vec [direc] :: direction

  • % wgt [real(4)] :: weight for thin sampling

  • % where [integer(2)] :: current obsSite no. (0 is initial value)

  • % asflag [integer(2)] :: non 0, if As has been generated from this ptcl (only for electrons)

  • % user [real(8)] :: user use

  • % inip [ptcl] :: Particle at production

  • % inipos [position]

  • % init [real(8)] :: time in length/beta (m)

  • % inivec [direc]

  • % parp [ptcl] :: parent particle

  • % parvec [direc]

  • % label [integer] :: (if LABELING > 0) put a label (1,2,…) on each particle. There is a global label_counter which is cleared at the start of 1 event generation. it is counted up when a particle is poped up from the stack. The label_counter is given to the label of the poped up particle. This may be needed to judge if the same particle crosses a given observation place more than once.

  • % info [integer] :: (if LABELING > 0) for each particle, when a particle is born this is initialized to 0. If the ptcl goes higher than 380km, 1 is added. This is for AMS observation.

Header :

Ztrack.h

type  element
Type fields:
  • % A [real(8)] :: mass number

  • % Z [real(8)] :: atomic number

  • % N [integer] :: nucleon number. N-Z is number of neutrons.

Header :

Zelement.h

type  epmedia
Type fields:
  • % noOfElem [integer] :: actual number of elements

  • % elem (maxElements) [element]

  • % No (maxElements) [real(8)] :: number of each element; copied to OrigNo and then normalized. so that sum be 1.0

  • % OrigNo (maxElements) [real(8)] :: number of each element

  • % MolMass [real(8)] :: olar mass normally no unit but in that case the value should be the same as the one in the unit of g/mol. default init value 0 will be given just before reading the data. Normally not used except for atmosphere. For the gas media, better to be given.

  • % w (maxElements) [real(8)] :: No(i)A(i)/sum(No(i)A(i)) same note as No

  • % npercm3 (maxElements) [real(8)] :: # of i-th element /cm3=No/Ai*rho*wi same note as No

  • % nsigma (maxElements) [real(8)] :: No(i)s_i

  • % sumns [real(8)] :: sum of above

  • % ndpsigma (maxElements) [real(8)] :: No(i)dps_i for DP

  • % sumndps [real(8)] :: sum of above

  • % sumNo [real(8)] :: sum of OrigNo(:)

  • % colElem [integer] :: element # at which interaction took place

  • % colA [integer] :: A of such one:. int(elem(colElm)%A+0.5)

  • % colZ [integer] :: Z of such one: int(elem(colElm)%Z)

  • % colXs [real(8)] :: x-section for that target(mb)

  • % xs [real(8)] :: this is xs for the media.

  • % ndensity [real(8)] :: effective number density /cm^3

  • % wp [real(8)] :: plasma frequency x hbar (GeV)

  • % n [real(8)] :: refractive index

  • % nd [real(8)] :: number of ingredients / cm^3

  • % A [real(8)] :: sum No x Ai

  • % Z [real(8)] :: sum No x Zi

  • % Z2 [real(8)] :: sum No x Zi**2

  • % ZZ1 [real(8)] :: sum No x Zi(Zi+1) for electron; this * t /(gamma beta2)^2 = Xc^2 (t g/cm2)

  • % MoliereForXc2 [real(8)]

  • % MoliereExpb [real(8)] :: exp(b) = t x this (t in g/cm2)a for z=1 and beta =1. this = 6702 sum/A sum = Sum No x Zi^(1/3)(Zi+1)/(1+3.327(Ziz/137)^2)

  • % Z1_3rd [real(8)] :: <Z^1/3> not <Z>^(1/3)

  • % Z2_3rd [real(8)] :: <Z^2/3> not <Z>^(2/3)

  • % mbtoPgrm [real(8)] :: 10^-27 x N0/A. If multiplied to sigma in mb, we obtain probability / (g/cm2).

  • % mbtoPkgrm [real(8)] :: mbtoPgrm/10d0

  • % mbtoPcm [real(8)] :: rho x mbtoPgrm. If multiplied to sigma in mb, we obtaind probability / cm

  • % mbtoPX0 [real(8)] :: mbtoPgrm x X0g. If multiplied to sigma in mb, we obtain probability /radation length. next ones are used when we approximate a compound /molecule as an atom

  • % mbtoPgrm2 [real(8)]

  • % mbtoPcm2 [real(8)]

  • % mbtoPX02 [real(8)]

  • % Z2byAeff [real(8)] :: sum wi x Zi**2/Ai

  • % Z5byAeff [real(8)] :: sum wi x Zi**5/Ai

  • % Aeff [real(8)] :: sum wi x Ai

  • % Z2eff [real(8)] :: Z2byAeff x Aeff

  • % Zeff [real(8)] :: sqrt(Z2eff)

  • % Zeff3 [real(8)] :: Zeff**(1/3)

  • % LogZ [real(8)] :: log(Zeff)

  • % A2eff [real(8)] :: sum wi x Ai^2

  • % ZbyAeff [real(8)] :: sum wi x Zi/Ai

  • % I [real(8)] :: average ionization potential energy in GeV.

  • % rho [real(8)] :: density in g/cm^3

  • % X0 [real(8)] :: radiation length. in cm

  • % X0m [real(8)] :: radiation lenght in m.

  • % X0g [real(8)] :: radiation length. in g/cm^2

  • % X0kg [real(8)] :: reaiation length in kg/m^2 =X0g*10

  • % gtocm [real(8)] :: g/cm^2 to cm.

  • % kgtom [real(8)] :: gtocm*1.0d-3: kg/m2 to m.

  • % dEdxatp3m [real(8)] :: dE/dx at p=3me for electron. ~ Ecrit

  • % Ecrit [real(8)] :: GeV for electron

  • % Ecritmu [real(8)] :: GeV for muon

  • % rhoc [real(8)] :: comp.rhoc is copied whenever new comp. comes note;this is real*8 while comp.rhoc is real*4

  • % gasF [integer] :: flag for gas. If 1, media is gas, 0 –>solid

  • % name [character(8)] :: name of media

  • % format [integer] :: format of the basic table. (1 or 2)

  • % s1 [real(8)] :: Migdal’s s1

  • % logs1 [real(8)] :: log(s1)

  • % basearea [real(8)] :: pi x Re**2 * N* Z /A *X0g = 0.15 Z/A*X0g

  • % cScrC1 [real(8)] :: const which appears in the complete screening crossection

  • % cScrC2 [real(8)] :: the other such one

  • % cScrMain [real(8)] :: (4/3C1 + C2)

  • % BirksC1 [real(8)] :: quenching correction coef.

  • % BirksC2 [real(8)]

  • % BirksCC [real(8)]

  • % Birks [character(1)] :: flag to identify what quenching correction should be applied using BirksC1, etc.

  • % srim [integer] :: index for srim data in module srimdata

  • % tbl [bpTbl]

  • % sh [sternh]

  • % cnst [SmpCnst]

  • % pe [photoE]

  • % urb [urban]

  • % mu [mubpn]

  • % xcom [epxcom]

Header :

Zmedia.h

type  SmpCnst
Type fields:
  • % CompScrE [real(8)] :: Energy above which we can use complete screening cross-sections. evaluate at Eg/Ee=x= 0.99.

  • % BrScrE [real(8)] :: below this, partial screened cross-section is needed ( = ComScrE)

  • % BremEgmin [real(8)] :: min. ratio Eg/Ee for brems at high energy region

  • % BremEemin [real(8)] :: Below this, partial screening brems x-section is not made. (Seltzer table 1 is used)

  • % BremLEemin [real(8)] :: log10 of BremEemin

  • % BremEeminLPM [real(8)] :: Min. energy of e+/e- above which LPM can be applied, if wanted.( Acutal application will be done if, Ee > Flpm*this and LPMeffect=T. Flpm and LPMeffect can be controled by epicsfile

  • % BremTXTL [integer] :: Size of the Brems total x-section table. in the energy region BremEemin ~ BrScrE: ~log10(BrScrE/BremEemin)*10

  • % BremEsize [integer] :: Size of log10 energy for 2D table for brems in the region A. BremTXTL/2

  • % BremUminLA [real(8)] :: min of uniform random number in the region A at energies BremEemin ~ BrScrE: 0.1

  • % BremUmaxLA [real(8)] :: max of uniform random number in the region A at energies BremEemin ~ BrScrE: 1.0

  • % BremUszLA [integer] :: Size of uniform random nubmbers for 2D table for brems in the region A.: 20

  • % BremdULA [real(8)] :: step of u in region A at Low energies

  • % BremdETXL [real(8)] :: log10E step for brem total cross secton log10(BrScrE/BremEemin)/(BremTXTL-1)

  • % BremdEL [real(8)] :: log10 step for 2D brem table at low energies

  • % BremUminLB [real(8)] :: min of uniform random number in the region B0. sqrt(u)

  • % BremUmaxLB [real(8)] :: max. sqrt(BremUminLA)

  • % BremUszLB [integer] :: u talbe size in region B. 20

  • % BremdULB [real(8)] :: step of u in B

  • % PairEgmin [real(8)] :: min. Eg above which pair cross section is computed 1.1 MeV. However, at energies from PairEgmin to PairNonSc, B.H original xsection is used as xsec= Norm * B.H where Norm * B.H(10MeV) = Pair(10MeV)

  • % PairNonSc [real(8)] :: see above.

  • % PairLEgmin [real(8)] :: log10 of PairEgmin

  • % PairEgmaxL [real(8)] :: Eg where LPM effect starts to appear

  • % PrScrE [real(8)] :: below this, screened cross-section is used.

  • % PairTXTL [integer] :: Size of the Pair total x-section table. in the energy region PairEgmin ~ PairEgmaxL; log10(PairEgmaxL/PairEgemin)*10

  • % PairEsize [integer] :: Size of log10 energy for 2D table for pair in the region A,B. PairTXTL/2

  • % PairUminLA [real(8)] :: min of uniform random number in the region A at energies PairEgmin ~ PairEgmaxL: 0.05

  • % PairUmaxLA [real(8)] :: max of uniform random number in the region A at energies PairEgmin ~ PairEgmaxL: 1.0

  • % PairUszLA [integer] :: Size of uniform random nubmbers for 2D table for pair in the region A.: 20

  • % PairdULA [real(8)] :: step of u in region A at Low energies

  • % PairdETXL [real(8)] :: log10 step of total pair cross-sec at low E

  • % PairUminLB [real(8)] :: min of uniform random number in the region B 0. sqrt(u)

  • % PairUmaxLB [real(8)] :: max. sqrt(PairUminLA)

  • % PairUszLB [integer] :: u talbe size in region B. 20

  • % PairdULB [real(8)] :: PairdULB=(PairUmaxLB-PairUminLB)/(PairUszLB-1))

  • % PairdELA [real(8)] :: log10(PairEgmaxL/ PairEgmin) /( PairEsize-1)

  • % PairdELB [real(8)] :: sqrt( log10(PairEgmaxL/ PairEgmin) )/( PairEsize-1)

  • % BrEeminS [real(8)] :: for Seltzer cross-section; lower energy region

  • % BrEgminS [real(8)] :: Eg min for Seltzer Brems. (not ratio 1keV)

  • % BrLEeminS [real(8)]

  • % BrEemaxS [real(8)]

  • % BrTXTS [integer]

  • % BrES [integer]

  • % BrUminSA [real(8)]

  • % BrUmaxSA [real(8)]

  • % BrUszSA [integer]

  • % BrdUSA [real(8)]

  • % BrdETXS [real(8)]

  • % BrdES [real(8)]

  • % BrUszSB [integer]

  • % BrUminSB [real(8)]

  • % BrUmaxSB [real(8)]

  • % BrdUSB [real(8)]

  • % BrEeminS2 [real(8)] :: for Seltzer cross-section; higher energy region upto 10 GeV

  • % BrEgminS2 [real(8)] :: Eg/Ee min for Seltzer Brems.

  • % BrLEeminS2 [real(8)]

  • % BrEemaxS2 [real(8)]

  • % BrTXTS2 [integer]

  • % BrES2 [integer]

  • % BrUminSA2 [real(8)]

  • % BrUmaxSA2 [real(8)]

  • % BrUszSA2 [integer]

  • % BrdUSA2 [real(8)]

  • % BrdETXS2 [real(8)]

  • % BrdES2 [real(8)]

  • % BrUszSB2 [integer]

  • % BrUminSB2 [real(8)]

  • % BrUmaxSB2 [real(8)]

  • % BrdUSB2 [real(8)]

  • % BrEgminH [real(8)] :: for LPM

  • % BrEe1H [real(8)]

  • % BrLEe1H [real(8)] :: log10( BrEe1H)

  • % BrneH [integer]

  • % BrdU1H [real(8)]

  • % BrdEH [real(8)] :: log E step cnst.BrdEH= log10(cnst.BrEe2H/cnst.BrEe1H)/(cnst.BrneH-1) inverse of the above

  • % BrEe2H [real(8)] :: max Ee where table is available

  • % BrU1H [real(8)]

  • % BrU2H [real(8)]

  • % Brnu1H [integer] :: =(cnst.BrU2H-cnst.BrU1H+0.00001d0)/cnst.BrdU1H+1

  • % BrneH2 [integer] :: for 2D table E size

  • % BrdEH2 [real(8)] :: // E bin

  • % BrEe2H2 [real(8)] :: max E for 2D table

  • % BrU3H [real(8)]

  • % BrU4H [real(8)]

  • % Brnu2H [integer]

  • % BrdVU2H [integer] :: = cnst.Brnu2H-1

  • % BrdU2H [real(8)] :: = (cnst.BrU4H - cnst.BrU3H)/cnst.BrdVU2H

  • % BrPow [real(8)]

  • % PrEg1H [real(8)] :: minimum Eg above which LPM works

  • % PrLEg1H [real(8)] :: log10 of PrEg1H

  • % PrneH [integer] :: number of Eg bins

  • % PrdU1H [real(8)] :: du

  • % PrdEH [real(8)] :: dE in log10(Eg)

  • % PrU1H [real(8)] :: minimum u= 0

  • % PrU2H [real(8)] :: maximum u= 1

  • % Prnu1H [integer] :: numboer of u bins

  • % PrEg2H [real(8)] :: max Eg where table is available. ———–muons nuclear interaction

  • % muNVmin [real(8)] :: min of Eg(virtual)/Emu by muon nuc. int.

  • % muNdU [real(8)] :: du for sampling table

  • % muNTXT [integer] :: total xs, dEdx(v<vmin), dEdx(vall), tab size.

  • % muNEmin [real(8)] :: above this, muon nuc. int. is treatable

  • % muNLEmin [real(8)] :: log10 of muNEmin

  • % muNEmax [real(8)] :: above this, use some scaling(sampling)

  • % muNEmax1 [real(8)] :: max E of 1D table

  • % muNdETX [real(8)] :: log10 Energy step for total muon nuc. int prob.

  • % muNdE [real(8)] :: log10 Energy step for sampling table

  • % muNUsize [integer] :: sampling table size for u.

  • % muNEsize [integer] :: sampling table size for log10 E

  • % muNpwtx [real(8)] :: prob/X0 energy dependence; power. set after table for total prob. is read

  • % muNpwdEdx0 [real(8)] :: dEdx(v<vmin)/Emu enery dependence; power set after table is read

  • % muNpwdEdxt [real(8)] :: dEdXt(v<vmax)/Emu energy dependence, power set after table is read

  • % muBrVmin [real(8)] :: brems min of Eg/Emu. for muon Brems

  • % muBrdU [real(8)] :: du for sampling table

  • % muBrTXT [integer] :: total xs, dEdx(v<vmin), dEdx(vall), tab size.

  • % muBrEmin [real(8)] :: above this, muon brems is treatable

  • % muBrLEmin [real(8)] :: log10 of muBrEmin

  • % muBrEmax [real(8)] :: above this, use some scaling

  • % muBrEmax1 [real(8)] :: max E of 1D table

  • % muBrdETX [real(8)] :: log10 Energy step for total muon brems prob.

  • % muBrdE [real(8)] :: log10 Energy step for sampling table

  • % muBrUsize [integer] :: sampling table size for u.

  • % muBrEsize [integer] :: sampling table size for log10 E dependence can be neglected

  • % muPrVmin [real(8)] :: pair creation min of Eg(virtual)/Emu by muon pair cre.

  • % muPrdU [real(8)] :: du for sampling table

  • % muPrTXT [integer] :: total xs, dEdx(v<vmin), dEdx(vall), tab size.

  • % muPrEmin [real(8)] :: above this, muon pair creation is treatable

  • % muPrLEmin [real(8)] :: log10 of muPrEmin

  • % muPrEmax [real(8)] :: above this, use some scaling

  • % muPrEmax1 [real(8)] :: max E of 1D table

  • % muPrdETX [real(8)] :: log10 Energy step for total muon pair prob.

  • % muPrdE [real(8)] :: log10 Energy step for sampling table

  • % muPrUsize [integer] :: sampling table size for u.

  • % muPrEsize [integer] :: sampling table size for log10 E dependence can be neglected

  • % how [integer]

  • % NormS [real(8)] :: normalization const

  • % NormPS [real(8)] :: normalization const

  • % NormCS [real(8)] :: normalization const

  • % NormSH [real(8)] :: normalization const

Header :

ZbpSample.h

type  bpTbl
Type fields:
  • % BrTXL (mxBrTXL,2) [real(8)]

  • % BrSTLA (mxBrTblLA, 1) [real(8)]

  • % BrSTLB (mxBrTblLB, 1) [real(8)]

  • % BrTXH (mxBrTXH,2) [real(8)]

  • % BrSTHA (mxBrTblHA, 1) [real(8)]

  • % BrSTHB (mxBrTblHB, 1) [real(8)]

  • % PrTXL (mxPrTXL) [real(8)]

  • % PrSTLA (mxPrTblLA, 1) [real(8)]

  • % PrSTLB (mxPrTblLB, 1) [real(8)]

  • % PrTXH (mxPrTXH) [real(8)]

  • % PrSTH (mxPrTblH, 1) [real(8)]

  • % BrTXS (mxBrTXS,2) [real(8)]

  • % BrSTSA (mxBrTblSA, 1) [real(8)]

  • % BrSTSB (mxBrTblSB, 1) [real(8)]

  • % BrTXS2 (mxBrTXS2,2) [real(8)]

  • % BrSTSA2 (mxBrTblSA2, 1) [real(8)]

  • % BrSTSB2 (mxBrTblSB2, 1) [real(8)]

  • % MuNTX (mxMuNTX) [real(8)]

  • % MuNdEdx0 (mxMuNTX) [real(8)]

  • % MuNdEdxt (mxMuNTX) [real(8)]

  • % MuBrTX (mxMuBrTX) [real(8)]

  • % MuBrdEdx0 (mxMuBrTX) [real(8)]

  • % MuBrdEdxt (mxMuBrTX) [real(8)]

  • % MuPrTX (mxMuPrTX) [real(8)]

  • % MuPrdEdx0 (mxMuPrTX) [real(8)]

  • % MuPrdEdxt (mxMuPrTX) [real(8)]

  • % MuNTbl (mxMuNTbl, 1) [real(8)]

  • % MuBrTbl (mxMuBrTbl, 1) [real(8)]

  • % MuPrTbl (mxMuPrTbl, 1) [real(8)]

Header :

ZbpTbl.h

type  site
Type fields:
  • % pos [position]

  • % Txyz2det (3,3) [real(8)] :: xyz to detector system transform mat

  • % Tdet2xyz (3,3) [real(8)] :: inverse of above

  • % zpl [real(8)] :: z value in 1ry system

  • % mu [real(8)] :: Moliere unit

  • % minitime [real(8)]

Header :

Zobsv.h

type  assite
Type fields:
  • % pos [position]

  • % zpl [real(8)]

  • % mu [real(8)] :: Moliere Unit

  • % esize [real(8)] :: electron size

  • % age [real(8)] :: size weighted age

Header :

Zobsv.h

type  magfield
Type fields:
  • % x [real(8)] :: in earth_center coordinate

  • % y [real(8)]

  • % z [real(8)]

  • % sys [character(4)] :: which system. ‘xyz’, ‘ned’, ‘hva’

Header :

Zmagfield.h

type  primaries
Type fields:
  • % each (maxNoOfComps) [component]

  • % cummInteFlux (maxNoOfComps) [real(8)]

  • % no_of_comps [integer] :: how many diff. compositions

  • % NoOfSamplings [integer] :: total number of samplings including ones discarded by cutoff

  • % NoOfSampComp (maxNoOfComps, 2) [integer] :: 1 is for number of sampling including discarded ones due to cutoff. 2 is only for employed ones. after a sampling of a 1ry, the following is fixed.

  • % label [integer] :: sampled primary label

  • % sampled_e [real(8)] :: sampled energy(or momentum) as defined in etype. If this is in total energy in GeV, it is the same as particle%fm%e

  • % particle [ptcl]

Header :

Zprimary.h

type  component

1 component of 1ry

Type fields:
  • % label [integer] :: composition label number

  • % symb [character(16)] :: ‘P’, ‘gamma’ etc.

  • % eunit [character(3)] :: ‘GeV’ etc

  • % etype [character(4)] :: ‘KE/n’ etc

  • % diff_or_inte [character(1)] :: ‘d’ or ‘i’

  • % flatterer [real(8)] :: dI/dE*E**flatterer

  • % cut [real(8)] :: lower cut off.

  • % cut2 [real(8)] :: upper cut off.

  • % energy (maxSegments+1) [real(8)] :: segment left energy

  • % flux (maxSegments+1) [real(8)] :: input flux dI/dE * E**flatterer above; from input table directly below; made by subroutines

  • % code [integer] :: particle code

  • % subcode [integer]

  • % charge [integer]

  • % togev [real(8)]

  • % norm_inte (maxSegments+1) [real(8)] :: normalized integral flux > E at segment left value

  • % beta (maxSegments+1) [real(8)] :: dI/dE=(true flux)= const*E**(-beta)

  • % no_of_seg [integer] :: no of segments given

  • % inte_value [real(8)] :: integral flux from min. E

  • % emin [real(8)] :: min and max energy defined

Header :

Zprimary.h

4.2. User accessible subroutines and functions

4.2.1. Useful

subroutine  cqEventNo(num, cumnum)

inquire the event number

Parameters:
  • num [integer,out] :: number of events in the current run

  • cumnum [integer,out] :: cummulative number of events so far. may be used after the initialization of an event, then this gives the number for that event. It will differ from num if Cont=t is used.

subroutine  cqIniRn(ir)

inquire the initial random seed. You may call this at any moment after the event initialization has been ended.

Parameters:

ir (2) [integer,out] :: the initial seed of the random number generator for the event.

subroutine  cqNoOfPrim(no)

inquire the number of primaries sampled specified in the primary spectrum file (MeV/n etc). You may call this after the event initialization has been ended. Note this is not necessary the total energy of the incident.

Parameters:

no [integer,out] :: no. of sampled primaries so far.

subroutine  cqPrimE(pOrE)

inquire sampled primary energy or p or rigidity as it is

Parameters:

pOrE [real(8),out] :: sampled primary energy or p or rigidity

subroutine  cqPrimary(prm)

inquire all about current primaries (input primary spectrum information)

Parameters:

prm [primaries,out] :: primary information

subroutine  cqFirstID(depth)

inquire the first interaction depth of the incident particle (vertical depth). After uv6.30, the knock-on process by p, He, etc is not regarded as the first interaction point; only their nuclear interaction is picked up. For non-nucleus, all interaction types are considered. If you want to have more detailed control, you may use chookNEPInt, chookGInt and/or chookEInt.

Parameters:

depth [real(8),out] :: the first interaction depth

subroutine  cqFirstIPI(aTrack)

inquire the complete first interaction point info. as a track of the primary. The definition of the first interaction point is the same as for cqFirstID.

Parameters:

aTrack [track,out] :: if ptrack%pos%depth=0. no interaction has been occurred

subroutine  cqIncident(aTrack, angleAtObs)

inquire incident particle. You may call this at any moment after the event initialization has been ended. The program unit containing this call must have #include "Ztrack.h".

Parameters:
  • aTrack [track,out] :: incident particle track information. You will see the ingredient by looking at the chookObs routine.

  • angleAtObs [coord,out] :: to get direction cosined of the incident particle in the “detector system” of the deepest observation level. angleAtObs%r(1), angleAtObs%r(2) and angleAtObs%r(3) are the three components.

4.2.2. Manager

subroutine  cwriteParam(io, force)

write parameters on the error output

Parameters:
  • io [integer,in] :: utput logical dev. #. ErrorOut –> stderr

  • force [integer,in] :: if non zero, Hidden parameters are written. hidden ones are also written when Hidden=T

subroutine  cprintPrim(out)

print primary information

Parameters:

out [integer,in] :: output logical device #

subroutine  cprintPrim(out)

print observation information

Parameters:

out [integer,in] :: output logical device #

subroutine  ckf2cos(kf, code, subcode, chg)

kf code to cosmos code.

Parameters:
  • kf [integer,in]

  • code [integer,out]

  • subcode [integer,out]

  • chg [integer,out]

subroutine  ccos2kf(code, subcode, chg, kf)

cosmos code to kf code;

Parameters:
  • code [integer,in]

  • subcode [integer,in]

  • chg [integer,in]

  • kf [integer,out]

subroutine  modCodeConv/ccos2pdg(aPtcl, pdgcode)

convert from Cosmos particle structure to PDG code

Parameters:
  • aPtcl [ptcl,in] :: Cosmos particle structure

  • pdgcode [integer,out] :: PDG particle code

4.2.3. EM

subroutine  epResetEcrit(io, name, newV, oldV, icon)

reset Crittical energy of a given media with “name”

Parameters:
  • io [integer,in] ::

    output message device #

    0

    some message is put as standard Fortran error message

    6

    some message is put as sysout.

    >0

    assume logical device is open with that number

    <0

    no message is put, but see next

  • name [character(*),in] :: media name such as “Air” if media with “name” is not found eroor message is

  • newV [real(8),in] :: new critical energy (GeV) new value is set to media%Ecrit

  • oldV [real(8),out] :: E crit so far defined. (GeV)

  • icon [integer,out] :: 0 if ok. -1 if some error

4.2.4. Tracking

subroutine  cavedEdx(eno, age, dedx)

Average dedx (\(2.2 \times 10^{-2}\) is set for the test)

Parameters:
  • eno [real(8),in] :: electon size

  • age [real(8),in] :: age of the shower

  • dedx [real(8),out] :: average dedx as defined in GeV/(kg/m2)

subroutine  cgetNmu(eth, nmu)

compute muon numbers this is not yet made. tentatively nmu = 0 is given.

Parameters:
  • eth [real(8),in] :: Threshold energy of muons. (GeV)

  • nmu (NoOfASSites) [real(8),out] :: output. number of muons E>eth

subroutine  cxyz2det(ly, a, b)

convert coord value in the “xyz” system into “det” system.

Parameters:
  • ly [integer(2),in] :: level # of the observation depth. Its origin is used to convert particle coordinate (‘a’ in E-xyz) into the detector coordinate (‘b’). Detector origin is the crossing point of 1ry direction and spherical surface at a given depth (height). Z-axis is vertical, X-axis is XaxisFromSouth (~90 deg normally magnetic East at the BaseL. The x-y plane is tangential to the spherical surface at the origin. ly is normally MovedTrack%where (integer(2))

  • a [coord,in] :: coord in ‘xyz’

  • b [coord,out] :: coord in ‘det’

subroutine  cdet2xyz(ly, a, b)

convert coord value in the “det” system into “det” system. See the parameter description of cxyz2det()

subroutine  cxyz2detD(ly, a, b)

convert coord value in the “xyz” system into “det” system for Direction cos. See the parameter description of cxyz2det()

subroutine  cdet2xyzD(ly, a, b)

convert coord value in the “det” system into “det” system for Direction cos. See the parameter description of cxyz2det()

subroutine  cdet2prim(ly, a, b)

xyz to primary system coord. conversion

Parameters:
  • ly [integer(2),in] :: base level # (xyz2prim case)->BaseL)

  • a [coord,in] :: input. coord. in ‘xyz’

  • b [coord,out] :: output. transformed coord. in ‘prim’

subroutine  cxyz2primD(a, b)

cxyz2primD: xyz to primary system for Direction cos. See the parameter description of cxyz2prim()

subroutine  csetPos(location)

set position information for a given xyz

Parameters:

location [position,inout] :: coord part of location is input.

subroutine  ciniTracking(incident)

inquire first int. point info.

Parameters:

incident [track,in]

subroutine  cqFirstColMedia(A, Z, xs)

retrns first col. element and Xsecion on it

Parameters:
  • A [integer,out]

  • Z [integer,out]

  • xs [real(8),out] :: Xsection

4.2.5. Atmosphere

subroutine  cllh2eCent(llh, xyz)

convert llh coordinates to E-xyz

Parameters:
  • llh [coord,in] :: containing data in latitude, longitude, height.

  • xyz [coord,out] ::

    The coordinate system is such that the origin is at the center of the earth

    x-axis

    directed to (0, 0) latitude and longitude.

    y-axis

    directed to (0, 90) latitude and longitude.

    z-axis

    directed to the north pole.

    xyz.r(1)

    x coordinate value in m

    xyz.r(2)

    y

    xyz.r(3)

    z

    note: xyz can be the same as llh. time component is unchanged

function  cvh2temp(vh)

temperature in Kelvin of ‘site’

Parameters:

vh [real(8),in] :: height in meter

Return:

cv2temp [real(8)] :: temperature in Kelvin

function  cvh2den(vh)

density of air

Parameters:

vh [real(8),in] :: height in meter

Return:

cvh2den [real(8)] :: density in kg/m3

4.2.6. Geomag

subroutine  cgeomag(yearin, llh, h, icon)
Parameters:
  • yearin [real(8),in] :: such as 1990.5

  • llh [coord,in] :: position around the earth. in ‘llh’ form is better. if not ‘llh’ conversion is done here.

  • h [magfield,out] :: magnetic field is set in the form of ‘ned’ (north, east-down). The unit is T.

  • icon [integer,out] ::

    0

    o.k

    1

    too heigh location. result could be suspicious.

    2

    input parameter wrong….

subroutine  ctransMagTo(sys, pos, a, b)

transform magnetic field components in one coordinate sytem to another.

a.sys \ sys

‘xyz’

‘hva’

‘ned’

‘xyz’

o

o

o

‘hva’

o

o

o

‘ned’

o

o

o

Parameters:
  • sys [character*(*),in] :: ‘xyz’, ‘hva’ or ‘ned’. the target coordinate system where magnetic filed is represented.

  • pos [coord,in] :: position where mag is given

  • a [magfield,in]

  • b [magfield,out] :: transformed component, b.sys=sys

subroutine csetMagField(sys, b1, b2, b3, b)

set Calculated magnetic field to /magfield/ b

Parameters:
  • sys [character(3),in] :: which system. ‘xyz’, ‘hva’, ‘ned’ etc

  • b1 [real(8)] :: 3 components of mag. in the system ‘sys’

  • b [magfield,out] :: /magfield/

4.2.7. General

function  klena(cha)

actual length of character string. note: if no character dec. is given for string, klena=0 will result. e.g. a=’abc’, klena(a). but klena(‘abc’) is ok.

Parameters:

cha [character*(*),in] :: character string

Return:

klena [integer] :: length of the string

4.3. User accessible variables

NoOfSites [integer]

No of particle observation sites (Zobsv.h) NoOfSites can be set upto 50 for the default configuration, which is defined by MAX_NO_OF_SITES in Zmaxdef.h.

NoOfASSites [integer]

maximum index for observation depths. (Zobsv.h) NoOfASSites can be set upto 50 for the default configuration, which is defined by MAX_NO_OF_AS_SITES in Zmaxdef.h.

CompASNe (i) [real(8)]

component A.S size produced by the input electron. For depths where this value is 0, avoid doing something here. (i=1, NoOfASSites; index for observation depths) (Zobsv.h)

CompASAge (i) [real(8)]

age of component A.S produced by the input electron. If this value is 2.0, the A.S is assumed to be very old and the CompASNe(i) is 0. You should skip treating deeper depths. (i=1, NoOfASSites ; index for observation depths) (Zobsv.h)

ObsSites (i) [site]

i=1, NoOfSites (Zobsv.h)

ASObsSites (i) [assite]

i=1, NoOfASSites (Zobsv.h)

Media (i) [epmedia]

Media information of Media_no=i (Zmedia.h)

MediaNo [integer]

current Media number. (ZmediaLoft.h)

TrackBefMove [track]

track before moved (Ztrackv.h)

MovedTrack [track]

contain track moved (Ztrackv.h)

AngleAtObsCopy [coord]

direction cos at the deepest Obssite in ‘det’ system (Zincidentv.h)

DcAtObsXyz [coord]

transfroamtion of AngleAtObsCopy to “xyz” system (Zincidentv.h)