IAP GITLAB

Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • AirShowerPhysics/corsika
  • rulrich/corsika
  • AAAlvesJr/corsika
  • Andre/corsika
  • arrabito/corsika
  • Nikos/corsika
  • olheiser73/corsika
  • AirShowerPhysics/papers/corsika
  • pranav/corsika
9 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 22474 deletions
Processes/Sibyll/rndm_dbl.f deleted 100644 → 0
View file @ 3deb039c
C***********************************************************************
C
C interface to PHOJET double precision random number generator
C for SIBYLL \FR'14
C
C***********************************************************************
DOUBLE PRECISION FUNCTION S_RNDM(IDUMMY)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DUMMY = dble(IDUMMY)
S_RNDM= PHO_RNDM(DUMMY)
END
C***********************************************************************
C
C initialization routine for double precision random number generator
C calls PHO_RNDIN \FR'14
C
C***********************************************************************
SUBROUTINE RND_INI
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON /RNDMGAS/ ISET
ISET = 0
CALL PHO_RNDIN(12,34,56,78)
END
DOUBLE PRECISION FUNCTION GASDEV(Idum)
C***********************************************************************
C Gaussian deviation
C***********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
IMPLICIT INTEGER(I-N)
COMMON /RNDMGAS/ ISET
SAVE
DATA ISET/0/
gasdev=idum
IF (ISET.EQ.0) THEN
1 V1=2.D0*S_RNDM(0)-1.D0
V2=2.D0*S_RNDM(1)-1.D0
R=V1**2+V2**2
IF(R.GE.1.D0)GO TO 1
FAC=SQRT(-2.D0*LOG(R)/R)
GSET=V1*FAC
GASDEV=V2*FAC
ISET=1
ELSE
GASDEV=GSET
ISET=0
ENDIF
RETURN
END
C***********************************************************************
DOUBLE PRECISION FUNCTION PHO_RNDM(DUMMY)
C***********************************************************************
C
C random number generator
C
C initialization by call to PHO_RNDIN needed!
C
C the algorithm is taken from
C G.Marsaglia, A.Zaman: 'Toward a unversal random number generator'
C Florida State Univ. preprint FSU-SCRI-87-70
C
C implementation by K. Hahn (Dec. 88), changed to include possibility
C of saving / reading generator registers to / from file (R.E. 10/98)
C
C generator should not depend on the hardware (if a real has
C at least 24 significant bits in internal representation),
C the period is about 2**144,
C
C internal registers:
C U(97),C,CD,CM,I,J - seed values as initialized in PHO_RNDIN
C
C
C***********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
COMMON /PORAND/ U(97),C,CD,CM,I,J
100 CONTINUE
RNDMI = DUMMY
RNDMI = U(I)-U(J)
IF ( RNDMI.LT.0.D0 ) RNDMI = RNDMI+1.D0
U(I) = RNDMI
I = I-1
IF ( I.EQ.0 ) I = 97
J = J-1
IF ( J.EQ.0 ) J = 97
C = C-CD
IF ( C.LT.0.D0 ) C = C+CM
RNDMI = RNDMI-C
IF ( RNDMI.LT.0.D0 ) RNDMI = RNDMI+1.D0
IF((ABS(RNDMI).LT.0.D0).OR.(ABS(RNDMI-1.D0).LT.1.D-10)) GOTO 100
PHO_RNDM = RNDMI
END
CDECK ID>, PHO_RNDIN
SUBROUTINE PHO_RNDIN(NA1,NA2,NA3,NB1)
C***********************************************************************
C
C initialization of PHO_RNDM, has to be called before using PHO_RNDM
C
C input:
C NA1,NA2,NA3,NB1 - values for initializing the generator
C NA? must be in 1..178 and not all 1;
C 12,34,56 are the standard values
C NB1 must be in 1..168;
C 78 is the standard value
C
C***********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
COMMON /PORAND/ U(97),C,CD,CM,I,J
MA1 = NA1
MA2 = NA2
MA3 = NA3
MB1 = NB1
I = 97
J = 33
DO 20 II2 = 1,97
S = 0.D0
T = 0.5D0
DO 10 II1 = 1,24
MAT = MOD(MOD(MA1*MA2,179)*MA3,179)
MA1 = MA2
MA2 = MA3
MA3 = MAT
MB1 = MOD(53*MB1+1,169)
IF ( MOD(MB1*MAT,64).GE.32 ) S = S+T
T = 0.5D0*T
10 CONTINUE
U(II2) = S
20 CONTINUE
C = 362436.D0/16777216.D0
CD = 7654321.D0/16777216.D0
CM = 16777213.D0/16777216.D0
END
CDECK ID>, PHO_RNDSI
SUBROUTINE PHO_RNDSI(UIN,CIN,CDIN,CMIN,IIN,JIN)
C***********************************************************************
C
C updates internal random number generator registers using
C registers given as arguments
C
C***********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
DIMENSION UIN(97)
COMMON /PORAND/ U(97),C,CD,CM,I,J
DO 10 KKK = 1,97
U(KKK) = UIN(KKK)
10 CONTINUE
C = CIN
CD = CDIN
CM = CMIN
I = IIN
J = JIN
END
CDECK ID>, PHO_RNDSO
SUBROUTINE PHO_RNDSO(UOUT,COUT,CDOUT,CMOUT,IOUT,JOUT)
C***********************************************************************
C
C copies internal registers from randon number generator
C to arguments
C
C***********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
DIMENSION UOUT(97)
COMMON /PORAND/ U(97),C,CD,CM,I,J
DO 10 KKK = 1,97
UOUT(KKK) = U(KKK)
10 CONTINUE
COUT = C
CDOUT = CD
CMOUT = CM
IOUT = I
JOUT = J
END
CDECK ID>, PHO_RNDTE
SUBROUTINE PHO_RNDTE(IO)
C***********************************************************************
C
C test of random number generator PHO_RNDM
C
C input:
C IO defines output
C 0 output only if an error is detected
C 1 output independend on an error
C
C uses PHO_RNDSI and PHO_RNDSO to bring the random number generator
C to same status as it had before the test run
C
C***********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
C input/output channels
INTEGER LI,LO
COMMON /POINOU/ LI,LO
DIMENSION UU(97)
DIMENSION U(6),X(6),D(6)
DATA U / 6533892.D0 , 14220222.D0 , 7275067.D0 ,
& 6172232.D0 , 8354498.D0 , 10633180.D0 /
CALL PHO_RNDSO(UU,CC,CCD,CCM,II,JJ)
CALL PHO_RNDIN(12,34,56,78)
DO 10 II1 = 1,20000
XX = PHO_RNDM(SD)
10 CONTINUE
SD = 0.D0
DO 20 II2 = 1,6
X(II2) = 4096.D0*(4096.D0*PHO_RNDM(XX))
D(II2) = X(II2)-U(II2)
SD = SD+ABS(D(II2))
20 CONTINUE
CALL PHO_RNDSI(UU,CC,CCD,CCM,II,JJ)
IF ((IO.EQ.1).OR.(ABS(SD).GT.0.D-10)) THEN
WRITE(LO,50) (U(I),X(I),D(I),I=1,6)
ENDIF
50 FORMAT(/,' PHO_RNDTE: test of the random number generator:',/,
& ' expected value calculated value difference',/,
& 6(F17.1,F20.1,F15.3,/),
& ' generator has the same status as before calling PHO_RNDTE',/)
END
CDECK ID>, PHO_RNDST
SUBROUTINE PHO_RNDST(MODE,FILENA)
C***********************************************************************
C
C read / write random number generator status from / to file
C
C input: MODE 1 read registers from file
C 2 dump registers to file
C
C FILENA file name
C
C***********************************************************************
IMPLICIT NONE
SAVE
INTEGER MODE
CHARACTER*(*) FILENA
C input/output channels
INTEGER LI,LO
COMMON /POINOU/ LI,LO
DOUBLE PRECISION UU,CC,CCD,CCM
DIMENSION UU(97)
INTEGER I,II,JJ
CHARACTER*80 CH_DUMMY
IF(MODE.EQ.1) THEN
WRITE(LO,'(/,1X,2A,A,/)') 'PHO_RNDST: ',
& 'reading random number registers from file ',FILENA
OPEN(12,FILE=FILENA,ERR=1010,STATUS='OLD')
READ(12,*,ERR=1010) CH_DUMMY
DO I=1,97
READ(12,*,ERR=1010) UU(I)
ENDDO
READ(12,*,ERR=1010) CC
READ(12,*,ERR=1010) CCD
READ(12,*,ERR=1010) CCM
READ(12,*,ERR=1010) II,JJ
CLOSE(12)
CALL PHO_RNDSI(UU,CC,CCD,CCM,II,JJ)
ELSE IF(MODE.EQ.2) THEN
WRITE(LO,'(/,1X,2A,A,/)') 'PHO_RNDST: ',
& 'dumping random number registers to file ',FILENA
OPEN(12,FILE=FILENA,ERR=1010,STATUS='UNKNOWN')
CALL PHO_RNDSO(UU,CC,CCD,CCM,II,JJ)
WRITE(12,'(1X,A)',ERR=1020) 'random number status registers:'
DO I=1,97
WRITE(12,'(1X,1P,E28.20)',ERR=1020) UU(I)
ENDDO
WRITE(12,'(1X,1P,E28.20)',ERR=1020) CC
WRITE(12,'(1X,1P,E28.20)',ERR=1020) CCD
WRITE(12,'(1X,1P,E28.20)',ERR=1020) CCM
WRITE(12,'(1X,2I4)',ERR=1020) II,JJ
CLOSE(12)
ELSE
WRITE(LO,'(/,1X,2A,I6,/)') 'PHO_RNDST: ',
& 'called with invalid mode, nothing done (mode)',MODE
ENDIF
RETURN
1010 CONTINUE
WRITE(LO,'(1X,2A,A,/)') 'PHO_RNDST: ',
& 'cannot open or read file ',FILENA
RETURN
1020 CONTINUE
WRITE(LO,'(1X,A,A,/)') 'PHO_RNDST: ',
& 'cannot open or write file ',FILENA
RETURN
END
C----------------------------------------
C standard generator
C----------------------------------------
REAL FUNCTION S_RNDM_std(IDUMMY)
C...Generator from the LUND montecarlo
C...Purpose: to generate random numbers uniformly distributed between
C...0 and 1, excluding the endpoints.
COMMON/LUDATR/MRLU(6),RRLU(100)
SAVE /LUDATR/
EQUIVALENCE (MRLU1,MRLU(1)),(MRLU2,MRLU(2)),(MRLU3,MRLU(3)),
&(MRLU4,MRLU(4)),(MRLU5,MRLU(5)),(MRLU6,MRLU(6)),
&(RRLU98,RRLU(98)),(RRLU99,RRLU(99)),(RRLU00,RRLU(100))
C... Initialize generation from given seed.
S_RNDM_std = real(idummy)
IF(MRLU2.EQ.0) THEN
IF (MRLU1 .EQ. 0) MRLU1 = 19780503 ! initial seed
IJ=MOD(MRLU1/30082,31329)
KL=MOD(MRLU1,30082)
I=MOD(IJ/177,177)+2
J=MOD(IJ,177)+2
K=MOD(KL/169,178)+1
L=MOD(KL,169)
DO 110 II=1,97
S=0.
T=0.5
DO 100 JJ=1,24
M=MOD(MOD(I*J,179)*K,179)
I=J
J=K
K=M
L=MOD(53*L+1,169)
IF(MOD(L*M,64).GE.32) S=S+T
T=0.5*T
100 CONTINUE
RRLU(II)=S
110 CONTINUE
TWOM24=1.
DO 120 I24=1,24
TWOM24=0.5*TWOM24
120 CONTINUE
RRLU98=362436.*TWOM24
RRLU99=7654321.*TWOM24
RRLU00=16777213.*TWOM24
MRLU2=1
MRLU3=0
MRLU4=97
MRLU5=33
ENDIF
C...Generate next random number.
130 RUNI=RRLU(MRLU4)-RRLU(MRLU5)
IF(RUNI.LT.0.) RUNI=RUNI+1.
RRLU(MRLU4)=RUNI
MRLU4=MRLU4-1
IF(MRLU4.EQ.0) MRLU4=97
MRLU5=MRLU5-1
IF(MRLU5.EQ.0) MRLU5=97
RRLU98=RRLU98-RRLU99
IF(RRLU98.LT.0.) RRLU98=RRLU98+RRLU00
RUNI=RUNI-RRLU98
IF(RUNI.LT.0.) RUNI=RUNI+1.
IF(RUNI.LE.0.OR.RUNI.GE.1.) GOTO 130
C...Update counters. Random number to output.
MRLU3=MRLU3+1
IF(MRLU3.EQ.1000000000) THEN
MRLU2=MRLU2+1
MRLU3=0
ENDIF
S_RNDM_std=RUNI
END
Processes/Sibyll/sibyll2.3c.cc deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/process/sibyll/sibyll2.3c.h>
#include <corsika/random/RNGManager.h>
#include <random>
int get_nwounded() { return s_chist_.nwd; }
double get_sibyll_mass2(int& id) { return s_mass1_.am2[id]; }
double s_rndm_(int&) {
static corsika::random::RNG& rng =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
std::uniform_real_distribution<double> dist;
return dist(rng);
}
Processes/Sibyll/sibyll2.3c.f deleted 100644 → 0
View file @ 3deb039c
Source diff could not be displayed: it is too large. Options to address this: view the blob.
Processes/Sibyll/sibyll2.3c.h deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_sib23c_interface_h_
#define _include_sib23c_interface_h_
//----------------------------------------------
// C++ interface for the SIBYLL event generator
//----------------------------------------------
// wrapper
extern "C" {
typedef char s_name[6];
// SIBYLL particle stack (FORTRAN COMMON)
// variables are: np : numer of particles on stack
// p : 4momentum + mass of particles on stack
// llist : id of particles on stack
extern struct {
double p[5][8000];
int llist[8000];
int np;
} s_plist_;
// additional information about interactions.
// number of wounded nucleons, number of hard and soft scatterings etc.
extern struct { int nnsof[20], nnjet[20], jdif[20], nwd, njet, nsof; } s_chist_;
extern struct {
double cbr[223 + 16 + 12 + 8];
int kdec[1338 + 6 * (16 + 12 + 8)];
int lbarp[99];
int idb[99];
} s_csydec_;
// additional particle stack for the mother particles of unstable particles
// stable particles have entry zero
extern struct { int llist1[8000]; } s_plist1_;
// tables with particle properties
// charge, strangeness and baryon number
extern struct {
int ichp[99];
int istr[99];
int ibar[99];
} s_chp_;
// tables with particle properties
// mass and mass squared
extern struct {
double am[99];
double am2[99];
} s_mass1_;
// table with particle names
extern struct { char namp[6][99]; } s_cnam_;
// debug info
extern struct {
int ncall;
int ndebug;
int lun;
} s_debug_;
// lund random generator setup
// extern struct {int mrlu[6]; float rrlu[100]; }ludatr_;
// sibyll main subroutine
void sibyll_(const int&, const int&, const double&);
// subroutine to initiate sibyll
void sibyll_ini_();
// subroutine to SET DECAYS
void dec_ini_();
// subroutine to initiate random number generator
// void rnd_ini_();
// print event
void sib_list_(int&);
// decay routine
void decsib_();
// interaction length
// double fpni_(double&, int&);
void sib_sigma_hnuc_(const int&, const int&, const double&, double&, double&, double&);
void sib_sigma_hp_(const int&, const double&, double&, double&, double&, double*, double&,
double&);
double s_rndm_(int&);
int get_nwounded();
double get_sibyll_mass2(int&);
// phojet random generator setup
void pho_rndin_(int&, int&, int&, int&);
}
#endif
Processes/Sibyll/sibyll_codes.dat deleted 100644 → 0
View file @ 3deb039c
# input file for particle conversion to/from SIBYLL
# the format of this file is: "corsika-identifier" "sibyll-id" "can-interact-in-sibyll" "cross-section-type"
Electron 3 0 0
Positron 2 0 0
NuE 15 0 0
NuEBar 16 0 0
MuMinus 5 0 0
MuPlus 4 0 0
NuMu 17 0 0
NuMuBar 18 0 0
TauMinus 91 0 0
TauPlus 90 0 0
NuTau 92 0 0
NuTauBar 93 0 0
Gamma 1 0 0
Pi0 6 1 2
# rho0 could interact but sibyll has no cross section/interaction length. was used for gamma had int
Rho0 27 0 0
K0Long 11 1 3
K0 21 0 3
K0Bar 22 0 3
PiPlus 7 1 2
PiMinus 8 1 2
RhoPlus 25 0 0
RhoMinus 26 0 0
Eta 23 0 0
Omega 32 0 0
K0Short 12 1 3
KStar0 30 0 0
KStar0Bar 31 0 0
KPlus 9 1 3
KMinus 10 1 3
KStarPlus 28 0 0
KStarMinus 29 0 0
DPlus 59 1 3
DMinus 60 1 3
DStarPlus 78 0 0
DStarMinus 79 0 0
D0 71 1 3
D0Bar 72 1 3
DStar0 80 0 0
DStar0Bar 81 0 0
DsPlus 74 1 3
DsMinus 75 1 3
DStarSPlus 76 0 0
DStarSMinus 77 0 0
EtaC 73 0 0
Neutron 14 1 1
AntiNeutron -14 1 1
Delta0 42 0 0
Delta0Bar -42 0 0
Proton 13 1 1
AntiProton -13 1 1
DeltaPlus 41 0 0
DeltaMinusBar -41 0 0
DeltaPlusPlus 40 0 0
DeltaMinusMinusBar -40 0 0
SigmaMinus 36 1 1
SigmaPlusBar -36 1 1
Lambda0 39 1 1
Lambda0Bar -39 1 1
Sigma0 35 1 1
Sigma0Bar -35 1 1
SigmaPlus 34 1 1
SigmaMinusBar -34 1 1
XiMinus 38 1 1
XiPlusBar -38 1 1
Xi0 37 1 1
Xi0Bar -37 1 1
OmegaMinus 49 0 0
OmegaPlusBar -49 0 0
SigmaC0 86 0 0
SigmaC0Bar -86 0 0
SigmaStarC0 96 0 0
SigmaStarC0Bar -96 0 0
LambdaCPlus 89 1 1
LambdaCMinusBar -89 1 1
XiC0 88 1 1
XiC0Bar -88 1 1
SigmaCPlus 85 0 0
SigmaCMinusBar -85 0 0
SigmaStarCPlus 95 0 0
SigmaStarCMinusBar -95 0 0
SigmaCPlusPlus 84 0 0
SigmaCMinusMinusBar -84 0 0
SigmaStarCPlusPlus 94 0 0
SigmaStarCMinusMinusBar -94 0 0
XiCPlus 87 1 1
XiCMinusBar -87 1 1
OmegaC0 99 0 0
OmegaC0Bar -99 0 0
Jpsi 83 0 0
Phi 33 0 0
#Unknown 0 0 0
Processes/Sibyll/testSibyll.cc deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/random/RNGManager.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalUnits.h>
#include <catch2/catch.hpp>
using namespace corsika;
using namespace corsika::process::sibyll;
TEST_CASE("Sibyll", "[processes]") {
SECTION("Sibyll -> Corsika") {
REQUIRE(particles::Electron::GetCode() ==
process::sibyll::ConvertFromSibyll(process::sibyll::SibyllCode::Electron));
}
SECTION("Corsika -> Sibyll") {
REQUIRE(process::sibyll::ConvertToSibyll(particles::Electron::GetCode()) ==
process::sibyll::SibyllCode::Electron);
REQUIRE(process::sibyll::ConvertToSibyllRaw(particles::Proton::GetCode()) == 13);
}
SECTION("canInteractInSibyll") {
REQUIRE(process::sibyll::CanInteract(particles::Proton::GetCode()));
REQUIRE(process::sibyll::CanInteract(particles::Code::XiCPlus));
REQUIRE_FALSE(process::sibyll::CanInteract(particles::Electron::GetCode()));
REQUIRE_FALSE(process::sibyll::CanInteract(particles::SigmaC0::GetCode()));
REQUIRE_FALSE(process::sibyll::CanInteract(particles::Nucleus::GetCode()));
REQUIRE_FALSE(process::sibyll::CanInteract(particles::Helium::GetCode()));
}
SECTION("cross-section type") {
REQUIRE(process::sibyll::GetSibyllXSCode(particles::Code::Electron) == 0);
REQUIRE(process::sibyll::GetSibyllXSCode(particles::Code::K0Long) == 3);
REQUIRE(process::sibyll::GetSibyllXSCode(particles::Code::SigmaPlus) == 1);
REQUIRE(process::sibyll::GetSibyllXSCode(particles::Code::PiMinus) == 2);
}
}
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/process/sibyll/sibyll2.3c.h>
using namespace corsika::units::si;
using namespace corsika::units;
TEST_CASE("SibyllInterface", "[processes]") {
// setup environment, geometry
environment::Environment<environment::IMediumModel> env;
auto& universe = *(env.GetUniverse());
auto theMedium =
environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
auto const* nodePtr = theMedium.get();
universe.AddChild(std::move(theMedium));
const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
SECTION("InteractionInterface") {
setup::Stack stack;
const HEPEnergyType E0 = 100_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
auto particle = stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Proton, E0, plab, pos, 0_ns});
particle.SetNode(nodePtr);
corsika::stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
Interaction model;
model.Init();
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
}
SECTION("NuclearInteractionInterface") {
setup::Stack stack;
const HEPEnergyType E0 = 400_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
auto particle =
stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
particles::Code::Nucleus, E0, plab, pos, 0_ns, 4, 2});
particle.SetNode(nodePtr);
corsika::stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
Interaction hmodel;
NuclearInteraction model(hmodel, env);
model.Init();
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
}
SECTION("DecayInterface") {
setup::Stack stack;
const HEPEnergyType E0 = 10_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
auto particle = stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Lambda0, E0, plab, pos, 0_ns});
corsika::stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
Decay model;
model.Init();
/*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(projectile);
// run checks
[[maybe_unused]] const TimeType time = model.GetLifetime(particle);
}
SECTION("DecayConfiguration") {
Decay model;
const std::vector<particles::Code> particleTestList = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
particles::Code::Lambda0Bar, particles::Code::NuE, particles::Code::D0Bar};
for (auto& pCode : particleTestList) {
model.SetUnstable(pCode);
// get state of sibyll internal config
REQUIRE(0 <= s_csydec_.idb[abs(process::sibyll::ConvertToSibyllRaw(pCode)) - 1]);
model.SetStable(pCode);
// get state of sibyll internal config
REQUIRE(0 >= s_csydec_.idb[abs(process::sibyll::ConvertToSibyllRaw(pCode)) - 1]);
}
}
}
Processes/StackInspector/CMakeLists.txt deleted 100644 → 0
View file @ 3deb039c
set (
MODEL_SOURCES
StackInspector.cc
)
set (
MODEL_HEADERS
StackInspector.h
)
set (
MODEL_NAMESPACE
corsika/process/stack_inspector
)
add_library (ProcessStackInspector STATIC ${MODEL_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessStackInspector ${MODEL_NAMESPACE} ${MODEL_HEADERS})
set_target_properties (
ProcessStackInspector
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
# PUBLIC_HEADER "${MODEL_HEADERS}"
)
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
ProcessStackInspector
CORSIKAcascade
CORSIKAunits
CORSIKAgeometry
CORSIKAsetup
CORSIKAlogging
)
target_include_directories (
ProcessStackInspector
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (
TARGETS ProcessStackInspector
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
# PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
)
# --------------------
# code unit testing
CORSIKA_ADD_TEST (testStackInspector)
target_link_libraries (
testStackInspector
ProcessStackInspector
CORSIKAgeometry
CORSIKAunits
CORSIKAtesting
)
Processes/StackInspector/StackInspector.cc deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/logging/Logger.h>
#include <corsika/setup/SetupTrajectory.h>
#include <chrono>
#include <iomanip>
#include <iostream>
#include <limits>
using namespace std;
using namespace corsika;
using namespace corsika::particles;
using namespace corsika::units::si;
using namespace corsika::process::stack_inspector;
template <typename TStack>
StackInspector<TStack>::StackInspector(const int vNStep, const bool vReportStack,
const HEPEnergyType vE0)
: StackProcess<StackInspector<TStack>>(vNStep)
, fReportStack(vReportStack)
, fE0(vE0)
, fStartTime(std::chrono::system_clock::now()) {}
template <typename TStack>
StackInspector<TStack>::~StackInspector() {}
template <typename TStack>
process::EProcessReturn StackInspector<TStack>::DoStack(const TStack& vS) {
[[maybe_unused]] int i = 0;
HEPEnergyType Etot = 0_GeV;
for (const auto& iterP : vS) {
HEPEnergyType E = iterP.GetEnergy();
Etot += E;
if (fReportStack) {
geometry::CoordinateSystem& rootCS = geometry::RootCoordinateSystem::GetInstance()
.GetRootCoordinateSystem(); // for printout
auto pos = iterP.GetPosition().GetCoordinates(rootCS);
cout << "StackInspector: i=" << setw(5) << fixed << (i++) << ", id=" << setw(30)
<< iterP.GetPID() << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, "
<< " pos=" << pos << " node = " << iterP.GetNode();
if (iterP.GetPID() == Code::Nucleus) cout << " nuc_ref=" << iterP.GetNucleusRef();
cout << endl;
}
}
auto const now = std::chrono::system_clock::now();
const std::chrono::duration<double> elapsed_seconds = now - fStartTime;
std::time_t const now_time = std::chrono::system_clock::to_time_t(now);
double const progress = (fE0 - Etot) / fE0;
double const eta_seconds = elapsed_seconds.count() / progress;
std::time_t const eta_time = std::chrono::system_clock::to_time_t(
fStartTime + std::chrono::seconds((int)eta_seconds));
cout << "StackInspector: "
<< " time=" << std::put_time(std::localtime(&now_time), "%T")
<< ", running=" << elapsed_seconds.count() << " seconds"
<< " (" << setw(3) << int(progress * 100) << "%)"
<< ", nStep=" << GetStep() << ", stackSize=" << vS.GetSize()
<< ", Estack=" << Etot / 1_GeV << " GeV"
<< ", ETA=" << std::put_time(std::localtime(&eta_time), "%T") << endl;
return process::EProcessReturn::eOk;
}
template <typename TStack>
void StackInspector<TStack>::Init() {
fReportStack = false;
fStartTime = std::chrono::system_clock::now();
}
#include <corsika/cascade/testCascade.h>
#include <corsika/setup/SetupStack.h>
template class process::stack_inspector::StackInspector<setup::Stack>;
template class process::stack_inspector::StackInspector<TestCascadeStack>;
Processes/StackInspector/StackInspector.h deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _Physics_StackInspector_StackInspector_h_
#define _Physics_StackInspector_StackInspector_h_
#include <corsika/process/StackProcess.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <chrono>
namespace corsika::process {
namespace stack_inspector {
template <typename TStack>
class StackInspector : public corsika::process::StackProcess<StackInspector<TStack>> {
typedef typename TStack::ParticleType Particle;
using corsika::process::StackProcess<StackInspector<TStack>>::GetStep;
public:
StackInspector(const int vNStep, const bool vReportStack,
const corsika::units::si::HEPEnergyType vE0);
~StackInspector();
void Init();
EProcessReturn DoStack(const TStack&);
/**
* To set a new E0, for example when a new shower event is started
*/
void SetE0(const corsika::units::si::HEPEnergyType vE0) { fE0 = vE0; }
private:
bool fReportStack;
corsika::units::si::HEPEnergyType fE0;
decltype(std::chrono::system_clock::now()) fStartTime;
};
} // namespace stack_inspector
} // namespace corsika::process
#endif
Processes/StackInspector/testStackInspector.cc deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <catch2/catch.hpp>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/cascade/testCascade.h>
using namespace corsika::units::si;
using namespace corsika::process::stack_inspector;
using namespace corsika;
using namespace corsika::geometry;
TEST_CASE("StackInspector", "[processes]") {
auto const& rootCS =
geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
geometry::Point const origin(rootCS, {0_m, 0_m, 0_m});
geometry::Vector<units::si::SpeedType::dimension_type> v(rootCS, 0_m / second,
0_m / second, 1_m / second);
geometry::Line line(origin, v);
geometry::Trajectory<geometry::Line> track(line, 10_s);
TestCascadeStack stack;
stack.Clear();
HEPEnergyType E0 = 100_GeV;
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Electron, E0,
corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
Point(rootCS, {0_m, 0_m, 10_km}), 0_ns});
SECTION("interface") {
StackInspector<TestCascadeStack> model(1, true, E0);
model.Init();
[[maybe_unused]] const process::EProcessReturn ret = model.DoStack(stack);
}
}
Processes/SwitchProcess/CMakeLists.txt deleted 100644 → 0
View file @ 3deb039c
set (
MODEL_HEADERS
SwitchProcess.h
)
set (
MODEL_NAMESPACE
corsika/process/switch_process
)
add_library (ProcessSwitch INTERFACE)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessSwitch ${MODEL_NAMESPACE} ${MODEL_HEADERS})
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
ProcessSwitch
INTERFACE
CORSIKAunits
CORSIKAprocesssequence
)
target_include_directories (
ProcessSwitch
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
# --------------------
# code unit testing
CORSIKA_ADD_TEST(testSwitchProcess)
target_link_libraries (
testSwitchProcess
ProcessSwitch
CORSIKAstackinterface
CORSIKAtesting
)
Processes/SwitchProcess/SwitchProcess.h deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _corsika_SwitchProcess_h
#define _corsika_SwitchProcess_h
#include <corsika/process/InteractionProcess.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process::switch_process {
/**
* This process provides an energy-based switch between two interaction processes P1 and
* P1. For energies below the threshold, P1 is invoked, otherwise P2. Both can be either
* single interaction processes or multiple ones combined in a ProcessSequence. A
* SwitchProcess itself will always be regarded as a distinct case when assembled into a
* (greater) ProcessSequence.
*/
template <class TLowEProcess, class THighEProcess>
class SwitchProcess : public BaseProcess<SwitchProcess<TLowEProcess, THighEProcess>> {
TLowEProcess& fLowEProcess;
THighEProcess& fHighEProcess;
units::si::HEPEnergyType const fThresholdEnergy;
public:
SwitchProcess(TLowEProcess& vLowEProcess, THighEProcess& vHighEProcess,
units::si::HEPEnergyType vThresholdEnergy)
: fLowEProcess(vLowEProcess)
, fHighEProcess(vHighEProcess)
, fThresholdEnergy(vThresholdEnergy) {}
void Init() {
fLowEProcess.Init();
fHighEProcess.Init();
}
template <typename TParticle>
corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& p) {
return 1 / GetInteractionLength(p);
}
template <typename TParticle>
units::si::GrammageType GetInteractionLength(TParticle& vParticle) {
if (vParticle.GetEnergy() < fThresholdEnergy) {
if constexpr (is_process_sequence_v<TLowEProcess>) {
return fLowEProcess.GetTotalInteractionLength(vParticle);
} else {
return fLowEProcess.GetInteractionLength(vParticle);
}
} else {
if constexpr (is_process_sequence_v<THighEProcess>) {
return fHighEProcess.GetTotalInteractionLength(vParticle);
} else {
return fHighEProcess.GetInteractionLength(vParticle);
}
}
}
// required to conform to ProcessSequence interface. We cannot just
// implement DoInteraction() because we want to call SelectInteraction
// in case a member process is a ProcessSequence.
template <typename TParticle, typename TSecondaries>
EProcessReturn SelectInteraction(
TParticle& vP, TSecondaries& vS,
[[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select,
corsika::units::si::InverseGrammageType& lambda_inv_count) {
if (vP.GetEnergy() < fThresholdEnergy) {
if constexpr (is_process_sequence_v<TLowEProcess>) {
return fLowEProcess.SelectInteraction(vP, vS, lambda_select, lambda_inv_count);
} else {
lambda_inv_count += fLowEProcess.GetInverseInteractionLength(vP);
// check if we should execute THIS process and then EXIT
if (lambda_select < lambda_inv_count) {
fLowEProcess.DoInteraction(vS);
return EProcessReturn::eInteracted;
} else {
return EProcessReturn::eOk;
}
}
} else {
if constexpr (is_process_sequence_v<THighEProcess>) {
return fHighEProcess.SelectInteraction(vP, vS, lambda_select, lambda_inv_count);
} else {
lambda_inv_count += fHighEProcess.GetInverseInteractionLength(vP);
// check if we should execute THIS process and then EXIT
if (lambda_select < lambda_inv_count) {
fHighEProcess.DoInteraction(vS);
return EProcessReturn::eInteracted;
} else {
return EProcessReturn::eOk;
}
}
}
}
};
} // namespace corsika::process::switch_process
#endif
Processes/SwitchProcess/testSwitchProcess.cc deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/process/switch_process/SwitchProcess.h>
#include <corsika/stack/SecondaryView.h>
#include <corsika/stack/Stack.h>
#include <corsika/units/PhysicalUnits.h>
#include <catch2/catch.hpp>
#include <algorithm>
#include <random>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units::si;
class TestStackData {
public:
// these functions are needed for the Stack interface
void Init() {}
void Clear() { fData.clear(); }
unsigned int GetSize() const { return fData.size(); }
unsigned int GetCapacity() const { return fData.size(); }
void Copy(int i1, int i2) { fData[i2] = fData[i1]; }
void Swap(int i1, int i2) { std::swap(fData[i1], fData[i2]); }
// custom data access function
void SetData(unsigned int i, HEPEnergyType v) { fData[i] = v; }
HEPEnergyType GetData(const unsigned int i) const { return fData[i]; }
// these functions are also needed by the Stack interface
void IncrementSize() { fData.resize(fData.size() + 1); }
void DecrementSize() {
if (fData.size() > 0) { fData.pop_back(); }
}
// custom private data section
private:
std::vector<HEPEnergyType> fData;
};
/**
* From static_cast of a StackIterator over entries in the
* TestStackData class you get and object of type
* TestParticleInterface defined here
*
* It provides Get/Set methods to read and write data to the "Data"
* storage of TestStackData obtained via
* "StackIteratorInterface::GetStackData()", given the index of the
* iterator "StackIteratorInterface::GetIndex()"
*
*/
template <typename StackIteratorInterface>
class TestParticleInterface
: public corsika::stack::ParticleBase<StackIteratorInterface> {
public:
using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData;
using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
/*
The SetParticleData methods are called for creating new entries
on the stack. You can specifiy various parametric versions to
perform this task:
*/
// default version for particle-creation from input data
void SetParticleData(const std::tuple<HEPEnergyType> v) { SetEnergy(std::get<0>(v)); }
void SetParticleData(TestParticleInterface<StackIteratorInterface>& /*parent*/,
std::tuple<HEPEnergyType> v) {
SetEnergy(std::get<0>(v));
}
// here are the fundamental methods for access to TestStackData data
void SetEnergy(HEPEnergyType v) { GetStackData().SetData(GetIndex(), v); }
HEPEnergyType GetEnergy() const { return GetStackData().GetData(GetIndex()); }
};
using SimpleStack = corsika::stack::Stack<TestStackData, TestParticleInterface>;
// see issue 161
#if defined(__clang__)
using StackTestView = corsika::stack::SecondaryView<TestStackData, TestParticleInterface>;
#elif defined(__GNUC__) || defined(__GNUG__)
using StackTestView = corsika::stack::MakeView<SimpleStack>::type;
#endif
auto constexpr kgMSq = 1_kg / (1_m * 1_m);
template <int N>
struct DummyProcess : InteractionProcess<DummyProcess<N>> {
void Init() {}
template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const {
return N * kgMSq;
}
template <typename TSecondaries>
corsika::process::EProcessReturn DoInteraction(TSecondaries& vSec) {
// to figure out which process was selected in the end, we produce N
// secondaries for DummyProcess<N>
for (int i = 0; i < N; ++i) {
vSec.AddSecondary(std::tuple<HEPEnergyType>{vSec.GetEnergy() / N});
}
return EProcessReturn::eOk;
}
};
using DummyLowEnergyProcess = DummyProcess<1>;
using DummyHighEnergyProcess = DummyProcess<2>;
using DummyAdditionalProcess = DummyProcess<3>;
TEST_CASE("SwitchProcess from InteractionProcess") {
DummyLowEnergyProcess low;
DummyHighEnergyProcess high;
DummyAdditionalProcess proc;
switch_process::SwitchProcess switchProcess(low, high, 1_TeV);
auto seq = switchProcess << proc;
SimpleStack stack;
SECTION("low energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV});
auto p = stack.GetNextParticle();
// low energy process returns 1 kg/m²
SECTION("interaction length") {
REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(1));
REQUIRE(seq.GetTotalInteractionLength(p) / kgMSq == Approx(3. / 4));
}
// low energy process creates 1 secondary
//~ SECTION("SelectInteraction") {
//~ typename SimpleStack::ParticleType theParticle =
//~ stack.GetNextParticle(); // as in corsika::Cascade
//~ StackTestView view(theParticle);
//~ auto projectile = view.GetProjectile();
//~ InverseGrammageType invLambda = 0 / kgMSq;
//~ switchProcess.SelectInteraction(p, projectile, 0.01 / kgMSq, invLambda);
//~ REQUIRE(view.GetSize() == 1);
//~ }
}
SECTION("high energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{4_TeV});
auto p = stack.GetNextParticle();
// high energy process returns 2 kg/m²
SECTION("interaction length") {
REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2));
REQUIRE(seq.GetTotalInteractionLength(p) / kgMSq == Approx(6. / 5));
}
// high energy process creates 2 secondaries
SECTION("SelectInteraction") {
typename SimpleStack::ParticleType theParticle =
stack.GetNextParticle(); // as in corsika::Cascade
StackTestView view(theParticle);
auto projectile = view.GetProjectile();
InverseGrammageType invLambda = 0 / kgMSq;
switchProcess.SelectInteraction(p, projectile, 0.01 / kgMSq, invLambda);
REQUIRE(view.GetSize() == 2);
}
}
}
TEST_CASE("SwitchProcess from ProcessSequence") {
DummyProcess<1> innerA;
DummyProcess<2> innerB;
DummyProcess<3> outer;
DummyProcess<4> additional;
auto seq = innerA << innerB;
switch_process::SwitchProcess switchProcess(seq, outer, 1_TeV);
auto completeSeq = switchProcess << additional;
SimpleStack stack;
SECTION("low energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV});
auto p = stack.GetNextParticle();
SECTION("interaction length") {
REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2. / 3));
REQUIRE(completeSeq.GetTotalInteractionLength(p) / kgMSq == Approx(4. / 7));
}
SECTION("SelectInteraction") {
std::vector<int> numberOfSecondaries;
for (int i = 0; i < 1000; ++i) {
typename SimpleStack::ParticleType theParticle =
stack.GetNextParticle(); // as in corsika::Cascade
StackTestView view(theParticle);
auto projectile = view.GetProjectile();
double r = i / 1000.;
InverseGrammageType invLambda = r * 7. / 4 / kgMSq;
InverseGrammageType accumulator = 0 / kgMSq;
completeSeq.SelectInteraction(p, projectile, invLambda, accumulator);
numberOfSecondaries.push_back(view.GetSize());
}
auto const mean =
std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) /
numberOfSecondaries.size();
REQUIRE(mean == Approx(12. / 7.).margin(0.01));
}
}
SECTION("high energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{3.0_TeV});
auto p = stack.GetNextParticle();
SECTION("interaction length") {
REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(3));
REQUIRE(completeSeq.GetTotalInteractionLength(p) / kgMSq == Approx(12. / 7.));
}
SECTION("SelectInteraction") {
std::vector<int> numberOfSecondaries;
for (int i = 0; i < 1000; ++i) {
typename SimpleStack::ParticleType theParticle =
stack.GetNextParticle(); // as in corsika::Cascade
StackTestView view(theParticle);
auto projectile = view.GetProjectile();
double r = i / 1000.;
InverseGrammageType invLambda = r * 7. / 12. / kgMSq;
InverseGrammageType accumulator = 0 / kgMSq;
completeSeq.SelectInteraction(p, projectile, invLambda, accumulator);
numberOfSecondaries.push_back(view.GetSize());
}
auto const mean =
std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) /
numberOfSecondaries.size();
REQUIRE(mean == Approx(24. / 7.).margin(0.01));
}
}
}
Processes/TrackWriter/CMakeLists.txt deleted 100644 → 0
View file @ 3deb039c
set (
MODEL_SOURCES
TrackWriter.cc
)
set (
MODEL_HEADERS
TrackWriter.h
)
set (
MODEL_NAMESPACE
corsika/process/track_writer
)
add_library (ProcessTrackWriter STATIC ${MODEL_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessTrackWriter ${MODEL_NAMESPACE} ${MODEL_HEADERS})
set_target_properties (
ProcessTrackWriter
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
# PUBLIC_HEADER "${MODEL_HEADERS}"
)
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
ProcessTrackWriter
CORSIKAunits
CORSIKAparticles
CORSIKAgeometry
CORSIKAsetup
)
target_include_directories (
ProcessTrackWriter
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (
TARGETS ProcessTrackWriter
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
# PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
)
# --------------------
# code unit testing
# CORSIKA_ADD_TEST(testNullModel)
#target_link_libraries (
# testNullModel ProcessNullModel
# CORSIKAsetup
# CORSIKAgeometry
# CORSIKAunits
# CORSIKAthirdparty # for catch2
# )
Processes/TrackWriter/TrackWriter.cc deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <limits>
using namespace corsika::setup;
using Particle = Stack::ParticleType;
using Track = Trajectory;
namespace corsika::process::track_writer {
void TrackWriter::Init() {
using namespace std::string_literals;
fFile.open(fFilename);
fFile << "# PID, E / eV, start coordinates / m, displacement vector to end / m "s
<< '\n';
}
template <>
process::EProcessReturn TrackWriter::DoContinuous(Particle& vP, Track& vT) {
using namespace units::si;
auto const start = vT.GetPosition(0).GetCoordinates();
auto const delta = vT.GetPosition(1).GetCoordinates() - start;
auto const pdg = static_cast<int>(particles::GetPDG(vP.GetPID()));
fFile << pdg << ' ' << vP.GetEnergy() / 1_eV << ' ' << start[0] / 1_m << ' '
<< start[1] / 1_m << ' ' << start[2] / 1_m << " " << delta[0] / 1_m << ' '
<< delta[1] / 1_m << ' ' << delta[2] / 1_m << '\n';
return process::EProcessReturn::eOk;
}
template <>
units::si::LengthType TrackWriter::MaxStepLength(Particle&, Track&) {
return units::si::meter * std::numeric_limits<double>::infinity();
}
} // namespace corsika::process::track_writer
Processes/TrackWriter/TrackWriter.h deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _Processes_track_writer_TrackWriter_h_
#define _Processes_track_writer_TrackWriter_h_
#include <corsika/process/ContinuousProcess.h>
#include <corsika/units/PhysicalUnits.h>
#include <fstream>
#include <string>
namespace corsika::process::track_writer {
class TrackWriter : public corsika::process::ContinuousProcess<TrackWriter> {
public:
TrackWriter(std::string const& filename)
: fFilename(filename) {}
void Init();
template <typename Particle, typename Track>
corsika::process::EProcessReturn DoContinuous(Particle&, Track&);
template <typename Particle, typename Track>
corsika::units::si::LengthType MaxStepLength(Particle&, Track&);
private:
std::string const fFilename;
std::ofstream fFile;
};
} // namespace corsika::process::track_writer
#endif
Processes/TrackingLine/CMakeLists.txt deleted 100644 → 0
View file @ 3deb039c
set (
MODEL_HEADERS
TrackingLine.h
)
set (
MODEL_SOURCES
TrackingLine.cc
)
set (
MODEL_NAMESPACE
corsika/process/tracking_line
)
add_library (ProcessTrackingLine STATIC ${MODEL_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessTrackingLine ${MODEL_NAMESPACE} ${MODEL_HEADERS})
set_target_properties (
ProcessTrackingLine
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
# PUBLIC_HEADER "${MODEL_HEADERS}"
)
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
ProcessTrackingLine
CORSIKAsetup
CORSIKAutilities
CORSIKAenvironment
CORSIKAunits
CORSIKAenvironment
CORSIKAgeometry
)
target_include_directories (
ProcessTrackingLine
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
# #-- -- -- -- -- -- -- -- -- --
# #code unit testing
CORSIKA_ADD_TEST (testTrackingLine)
target_link_libraries (
testTrackingLine
ProcessTrackingLine
CORSIKAtesting
)
Processes/TrackingLine/TrackingLine.cc deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/environment/Environment.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/QuantityVector.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Vector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <limits>
#include <stdexcept>
#include <utility>
using namespace corsika::geometry;
using namespace corsika::units::si;
namespace corsika::process::tracking_line {
std::optional<std::pair<TimeType, TimeType>> TimeOfIntersection(Line const& line,
Sphere const& sphere) {
auto const delta = line.GetR0() - sphere.GetCenter();
auto const v = line.GetV0();
auto const vSqNorm =
v.squaredNorm(); // todo: get rid of this by having V0 normalized always
auto const R = sphere.GetRadius();
auto const vDotDelta = v.dot(delta);
auto const discriminant =
vDotDelta * vDotDelta - vSqNorm * (delta.squaredNorm() - R * R);
if (discriminant.magnitude() > 0) {
auto const sqDisc = sqrt(discriminant);
auto const invDenom = 1 / vSqNorm;
return std::make_pair((-vDotDelta - sqDisc) * invDenom,
(-vDotDelta + sqDisc) * invDenom);
} else {
return {};
}
}
TimeType TimeOfIntersection(Line const& vLine, Plane const& vPlane) {
auto const delta = vPlane.GetCenter() - vLine.GetR0();
auto const v = vLine.GetV0();
auto const n = vPlane.GetNormal();
auto const c = n.dot(v);
if (c.magnitude() == 0) {
return std::numeric_limits<TimeType::value_type>::infinity() * 1_s;
} else {
return n.dot(delta) / c;
}
}
} // namespace corsika::process::tracking_line
Processes/TrackingLine/TrackingLine.h deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_corsika_processes_TrackingLine_h_
#define _include_corsika_processes_TrackingLine_h_
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Plane.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <optional>
#include <type_traits>
#include <utility>
namespace corsika::environment {
template <typename IEnvironmentModel>
class Environment;
template <typename IEnvironmentModel>
class VolumeTreeNode;
} // namespace corsika::environment
namespace corsika::process {
namespace tracking_line {
std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
TimeOfIntersection(geometry::Line const&, geometry::Sphere const&);
corsika::units::si::TimeType TimeOfIntersection(geometry::Line const&,
geometry::Plane const&);
class TrackingLine {
public:
TrackingLine(){};
template <typename Particle> // was Stack previously, and argument was
// Stack::StackIterator
auto GetTrack(Particle const& p) {
using namespace corsika::units::si;
using namespace corsika::geometry;
geometry::Vector<SpeedType::dimension_type> const velocity =
p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
auto const currentPosition = p.GetPosition();
std::cout << "TrackingLine pid: " << p.GetPID()
<< " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
std::cout << "TrackingLine pos: "
<< currentPosition.GetCoordinates()
// << " [" << p.GetNode()->GetModelProperties().GetName() << "]"
<< std::endl;
std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
<< " GeV " << std::endl;
std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
// to do: include effect of magnetic field
geometry::Line line(currentPosition, velocity);
auto const* currentLogicalVolumeNode = p.GetNode();
//~ auto const* currentNumericalVolumeNode =
//~ fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
auto const numericallyInside =
currentLogicalVolumeNode->GetVolume().Contains(currentPosition);
std::cout << "numericallyInside = " << (numericallyInside ? "true" : "false");
auto const& children = currentLogicalVolumeNode->GetChildNodes();
auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes();
std::vector<std::pair<TimeType, decltype(p.GetNode())>> intersections;
// for entering from outside
auto addIfIntersects = [&](auto const& vtn) {
auto const& volume = vtn.GetVolume();
auto const& sphere = dynamic_cast<geometry::Sphere const&>(
volume); // for the moment we are a bit bold here and assume
// everything is a sphere, crashes with exception if not
if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
auto const [t1, t2] = *opt;
std::cout << "intersection times: " << t1 / 1_s << "; "
<< t2 / 1_s
// << " " << vtn.GetModelProperties().GetName()
<< std::endl;
if (t1.magnitude() > 0)
intersections.emplace_back(t1, &vtn);
else if (t2.magnitude() > 0)
std::cout << "inside other volume" << std::endl;
}
};
for (auto const& child : children) { addIfIntersects(*child); }
for (auto const* ex : excluded) { addIfIntersects(*ex); }
{
auto const& sphere = dynamic_cast<geometry::Sphere const&>(
currentLogicalVolumeNode->GetVolume());
// for the moment we are a bit bold here and assume
// everything is a sphere, crashes with exception if not
[[maybe_unused]] auto const [t1, t2] = *TimeOfIntersection(line, sphere);
[[maybe_unused]] auto dummy_t1 = t1;
intersections.emplace_back(t2, currentLogicalVolumeNode->GetParent());
}
auto const minIter = std::min_element(
intersections.cbegin(), intersections.cend(),
[](auto const& a, auto const& b) { return a.first < b.first; });
TimeType min;
if (minIter == intersections.cend()) {
min = 1_s; // todo: do sth. more reasonable as soon as tracking is able
// to handle the numerics properly
throw std::runtime_error("no intersection with anything!");
} else {
min = minIter->first;
}
std::cout << " t-intersect: "
<< min
// << " " << minIter->second->GetModelProperties().GetName()
<< std::endl;
return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
velocity.norm() * min, minIter->second);
}
};
} // namespace tracking_line
} // namespace corsika::process
#endif
Processes/TrackingLine/testTrackingLine.cc deleted 100644 → 0
View file @ 3deb039c
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/process/tracking_line/TrackingLine.h>
#include <testTrackingLineStack.h> // test-build, and include file is obtained from CMAKE_CURRENT_SOURCE_DIR
#include <corsika/environment/Environment.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Vector.h>
#include <corsika/setup/SetupTrajectory.h>
using corsika::setup::Trajectory;
#include <catch2/catch.hpp>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::geometry;
#include <iostream>
using namespace std;
using namespace corsika::units::si;
TEST_CASE("TrackingLine") {
environment::Environment<environment::Empty> env; // dummy environment
auto const& cs = env.GetCoordinateSystem();
tracking_line::TrackingLine tracking;
SECTION("intersection with sphere") {
Point const origin(cs, {0_m, 0_m, -5_m});
Point const center(cs, {0_m, 0_m, 10_m});
Sphere const sphere(center, 1_m);
Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
0_m / second, 1_m / second);
Line line(origin, v);
setup::Trajectory traj(line, 12345_s);
auto const opt =
tracking_line::TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
REQUIRE(opt.has_value());
auto [t1, t2] = opt.value();
REQUIRE(t1 / 14_s == Approx(1));
REQUIRE(t2 / 16_s == Approx(1));
auto const optNoIntersection =
tracking_line::TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
REQUIRE_FALSE(optNoIntersection.has_value());
}
SECTION("maximally possible propagation") {
auto& universe = *(env.GetUniverse());
auto const radius = 20_m;
auto theMedium = environment::Environment<environment::Empty>::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
auto const* theMediumPtr = theMedium.get();
universe.AddChild(std::move(theMedium));
TestTrackingLineStack stack;
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::MuPlus,
1_GeV,
{cs, {0_GeV, 0_GeV, 1_GeV}},
{cs, {0_m, 0_m, 0_km}},
0_ns});
auto p = stack.GetNextParticle();
p.SetNode(theMediumPtr);
Point const origin(cs, {0_m, 0_m, 0_m});
Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
0_m / second, 1_m / second);
Line line(origin, v);
auto const [traj, geomMaxLength, nextVol] = tracking.GetTrack(p);
[[maybe_unused]] auto dummy_geomMaxLength = geomMaxLength;
[[maybe_unused]] auto dummy_nextVol = nextVol;
REQUIRE((traj.GetPosition(1.) - Point(cs, 0_m, 0_m, radius))
.GetComponents(cs)
.norm()
.magnitude() == Approx(0).margin(1e-4));
}
}

data privacy - impressum