IAP GITLAB
Skip to content
GitLab
Explore
Sign in
Register
Primary navigation
Search or go to…
Project
C
corsika
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Build
Pipelines
Jobs
Pipeline schedules
Artifacts
Deploy
Releases
Package Registry
Model registry
Operate
Environments
Terraform modules
Monitor
Incidents
Analyze
Value stream analytics
Contributor analytics
CI/CD analytics
Repository analytics
Issue analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
Community forum
Contribute to GitLab
Provide feedback
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
Pranav Sampathkumar
corsika
Commits
bc34c375
Commit
bc34c375
authored
5 years ago
by
ralfulrich
Committed by
Ralf Ulrich
5 years ago
Browse files
Options
Downloads
Patches
Plain Diff
clang format
parent
20b59cc4
No related branches found
Branches containing commit
No related tags found
Tags containing commit
No related merge requests found
Changes
2
Hide whitespace changes
Inline
Side-by-side
Showing
2 changed files
Processes/QGSJetII/Interaction.cc
+86
-89
86 additions, 89 deletions
Processes/QGSJetII/Interaction.cc
Processes/QGSJetII/Interaction.h
+4
-4
4 additions, 4 deletions
Processes/QGSJetII/Interaction.h
with
90 additions
and
93 deletions
Processes/QGSJetII/Interaction.cc
+
86
−
89
View file @
bc34c375
...
...
@@ -64,12 +64,10 @@ namespace corsika::process::qgsjetII {
}
}
units
::
si
::
CrossSectionType
Interaction
::
GetCrossSection
(
const
particles
::
Code
BeamId
,
const
particles
::
Code
TargetId
,
const
units
::
si
::
HEPEnergyType
Elab
,
const
unsigned
int
Abeam
,
const
unsigned
int
Atarget
)
const
{
units
::
si
::
CrossSectionType
Interaction
::
GetCrossSection
(
const
particles
::
Code
BeamId
,
const
particles
::
Code
TargetId
,
const
units
::
si
::
HEPEnergyType
Elab
,
const
unsigned
int
Abeam
,
const
unsigned
int
Atarget
)
const
{
using
namespace
units
::
si
;
double
sigProd
=
std
::
numeric_limits
<
double
>::
infinity
();
...
...
@@ -92,9 +90,9 @@ namespace corsika::process::qgsjetII {
throw
std
::
runtime_error
(
"QgsjetII target outside range. "
);
}
cout
<<
"QgsjetII::GetCrossSection Elab="
<<
Elab
<<
" iBeam="
<<
iBeam
<<
" iProjectile="
<<
iProjectile
<<
" iTarget="
<<
iTarget
<<
endl
;
sigProd
=
qgsect_
(
Elab
/
1
_GeV
,
iBeam
,
iProjectile
,
iTarget
);
cout
<<
"QgsjetII::GetCrossSection Elab="
<<
Elab
<<
" iBeam="
<<
iBeam
<<
"
iProjectile="
<<
iProjectile
<<
"
iTarget="
<<
iTarget
<<
endl
;
sigProd
=
qgsect_
(
Elab
/
1
_GeV
,
iBeam
,
iProjectile
,
iTarget
);
cout
<<
"QgsjetII::GetCrossSection sigProd="
<<
sigProd
<<
endl
;
}
...
...
@@ -298,8 +296,7 @@ namespace corsika::process::qgsjetII {
}
cout
<<
"Interaction: "
<<
" DoInteraction: E(GeV):"
<<
eProjectileLab
/
1
_GeV
<<
endl
;
<<
" DoInteraction: E(GeV):"
<<
eProjectileLab
/
1
_GeV
<<
endl
;
count_
++
;
qgini_
(
eProjectileLab
/
1
_GeV
,
kBeam
,
projQgsCode
,
targetQgsCode
);
// this is from CRMC, is this REALLY needed ???
...
...
@@ -313,97 +310,97 @@ namespace corsika::process::qgsjetII {
// fragments
QGSJetIIFragmentsStack
qfs
;
for
(
auto
&
fragm
:
qfs
)
{
particles
::
Code
idFragm
=
particles
::
Code
::
Nucleus
;
int
A
=
fragm
.
GetFragmentSize
();
int
Z
=
0
;
switch
(
A
)
{
case
1
:
{
// proton/neutron
idFragm
=
particles
::
Code
::
Proton
;
MomentumVector
Plab_nucleon
(
rootCS
,
{
0.0
_GeV
,
0.0
_GeV
,
sqrt
((
eProjectileLab
+
particles
::
Proton
::
GetMass
())
*
(
eProjectileLab
-
particles
::
Proton
::
GetMass
()))});
auto
const
PlabRot
=
boost
.
fromCoM
(
FourVector
(
eProjectileLab
,
Plab_nucleon
));
cout
<<
"secondary fragment>"
//<< static_cast<int>(psec.GetPID())
<<
" "
<<
idFragm
<<
" eProjectileLab="
<<
eProjectileLab
<<
" "
<<
endl
;
auto
pnew
=
vP
.
AddSecondary
(
tuple
<
particles
::
Code
,
units
::
si
::
HEPEnergyType
,
stack
::
MomentumVector
,
geometry
::
Point
,
units
::
si
::
TimeType
>
{
idFragm
,
PlabRot
.
GetTimeLikeComponent
(),
PlabRot
.
GetSpaceLikeComponents
(),
pOrig
,
tOrig
});
Plab_final
+=
pnew
.
GetMomentum
();
Elab_final
+=
pnew
.
GetEnergy
();
}
break
;
case
2
:
// deuterium
Z
=
1
;
break
;
case
3
:
// tritium
Z
=
1
;
break
;
case
4
:
// helium
Z
=
2
;
break
;
default
:
// nucleus
{
Z
=
int
(
A
/
2.15
+
0.7
);
}
}
if
(
idFragm
==
particles
::
Code
::
Nucleus
)
{
MomentumVector
Plab_nucleus
(
rootCS
,
{
0.0
_GeV
,
0.0
_GeV
,
sqrt
((
eProjectileLab
+
constants
::
nucleonMass
*
A
)
*
(
eProjectileLab
-
constants
::
nucleonMass
*
A
))});
auto
const
PlabRot
=
boost
.
fromCoM
(
FourVector
(
eProjectileLab
*
A
,
Plab_nucleus
));
cout
<<
"secondary fragment>"
<<
" "
<<
idFragm
<<
" eProjectileLab="
<<
eProjectileLab
<<
" A="
<<
A
<<
" Z="
<<
Z
<<
endl
;
particles
::
Code
idFragm
=
particles
::
Code
::
Nucleus
;
int
A
=
fragm
.
GetFragmentSize
();
int
Z
=
0
;
switch
(
A
)
{
case
1
:
{
// proton/neutron
idFragm
=
particles
::
Code
::
Proton
;
MomentumVector
Plab_nucleon
(
rootCS
,
{
0.0
_GeV
,
0.0
_GeV
,
sqrt
((
eProjectileLab
+
particles
::
Proton
::
GetMass
())
*
(
eProjectileLab
-
particles
::
Proton
::
GetMass
()))});
auto
const
PlabRot
=
boost
.
fromCoM
(
FourVector
(
eProjectileLab
,
Plab_nucleon
));
cout
<<
"secondary fragment>"
//<< static_cast<int>(psec.GetPID())
<<
" "
<<
idFragm
<<
" eProjectileLab="
<<
eProjectileLab
<<
" "
<<
endl
;
auto
pnew
=
vP
.
AddSecondary
(
tuple
<
particles
::
Code
,
units
::
si
::
HEPEnergyType
,
stack
::
MomentumVector
,
geometry
::
Point
,
units
::
si
::
TimeType
,
unsigned
short
,
unsigned
short
>
{
idFragm
,
PlabRot
.
GetTimeLikeComponent
(),
PlabRot
.
GetSpaceLikeComponents
(),
pOrig
,
tOrig
,
A
,
Z
});
geometry
::
Point
,
units
::
si
::
TimeType
>
{
idFragm
,
PlabRot
.
GetTimeLikeComponent
(),
PlabRot
.
GetSpaceLikeComponents
(),
pOrig
,
tOrig
});
Plab_final
+=
pnew
.
GetMomentum
();
Elab_final
+=
pnew
.
GetEnergy
();
}
break
;
case
2
:
// deuterium
Z
=
1
;
break
;
case
3
:
// tritium
Z
=
1
;
break
;
case
4
:
// helium
Z
=
2
;
break
;
default
:
// nucleus
{
Z
=
int
(
A
/
2.15
+
0.7
);
}
}
// secondaries
QGSJetIIStack
qs
;
for
(
auto
&
psec
:
qs
)
{
auto
const
Plab
=
psec
.
GetMomentum
();
auto
const
Elab
=
psec
.
GetEnergy
();
// transform energy to lab. frame
auto
const
PlabRot
=
boost
.
fromCoM
(
FourVector
(
Elab
,
Plab
));
cout
<<
"secondary hadron>"
//<< static_cast<int>(psec.GetPID())
<<
" "
<<
process
::
qgsjetII
::
ConvertFromQgsjetII
(
psec
.
GetPID
())
<<
" Elab="
<<
Elab
<<
" "
<<
endl
;
// add to corsika stack
if
(
idFragm
==
particles
::
Code
::
Nucleus
)
{
MomentumVector
Plab_nucleus
(
rootCS
,
{
0.0
_GeV
,
0.0
_GeV
,
sqrt
((
eProjectileLab
+
constants
::
nucleonMass
*
A
)
*
(
eProjectileLab
-
constants
::
nucleonMass
*
A
))});
auto
const
PlabRot
=
boost
.
fromCoM
(
FourVector
(
eProjectileLab
*
A
,
Plab_nucleus
));
cout
<<
"secondary fragment>"
<<
" "
<<
idFragm
<<
" eProjectileLab="
<<
eProjectileLab
<<
" A="
<<
A
<<
" Z="
<<
Z
<<
endl
;
auto
pnew
=
vP
.
AddSecondary
(
tuple
<
particles
::
Code
,
units
::
si
::
HEPEnergyType
,
stack
::
MomentumVector
,
geometry
::
Point
,
units
::
si
::
TimeType
>
{
process
::
qgsjetII
::
ConvertFromQgsjetII
(
psec
.
GetPID
()),
PlabRot
.
GetTimeLikeComponent
(),
PlabRot
.
GetSpaceLikeComponents
(),
pOrig
,
tOrig
});
geometry
::
Point
,
units
::
si
::
TimeType
,
unsigned
short
,
unsigned
short
>
{
idFragm
,
PlabRot
.
GetTimeLikeComponent
(),
PlabRot
.
GetSpaceLikeComponents
(),
pOrig
,
tOrig
,
A
,
Z
});
Plab_final
+=
pnew
.
GetMomentum
();
Elab_final
+=
pnew
.
GetEnergy
();
}
cout
<<
"conservation (all GeV): Ecm_final= n/a"
/* << Ecm_final / 1_GeV*/
<<
endl
<<
"Elab_final="
<<
Elab_final
/
1
_GeV
<<
", Plab_final="
<<
(
Plab_final
/
1
_GeV
).
GetComponents
()
<<
", N_wounded,targ="
<<
QGSJetIIFragmentsStackData
::
GetWoundedNucleonsTarget
()
<<
", N_wounded,proj="
<<
QGSJetIIFragmentsStackData
::
GetWoundedNucleonsProjectile
()
<<
", N_fragm,proj="
<<
qfs
.
GetSize
()
<<
endl
;
}
// secondaries
QGSJetIIStack
qs
;
for
(
auto
&
psec
:
qs
)
{
auto
const
Plab
=
psec
.
GetMomentum
();
auto
const
Elab
=
psec
.
GetEnergy
();
// transform energy to lab. frame
auto
const
PlabRot
=
boost
.
fromCoM
(
FourVector
(
Elab
,
Plab
));
cout
<<
"secondary hadron>"
//<< static_cast<int>(psec.GetPID())
<<
" "
<<
process
::
qgsjetII
::
ConvertFromQgsjetII
(
psec
.
GetPID
())
<<
" Elab="
<<
Elab
<<
" "
<<
endl
;
// add to corsika stack
auto
pnew
=
vP
.
AddSecondary
(
tuple
<
particles
::
Code
,
units
::
si
::
HEPEnergyType
,
stack
::
MomentumVector
,
geometry
::
Point
,
units
::
si
::
TimeType
>
{
process
::
qgsjetII
::
ConvertFromQgsjetII
(
psec
.
GetPID
()),
PlabRot
.
GetTimeLikeComponent
(),
PlabRot
.
GetSpaceLikeComponents
(),
pOrig
,
tOrig
});
Plab_final
+=
pnew
.
GetMomentum
();
Elab_final
+=
pnew
.
GetEnergy
();
}
cout
<<
"conservation (all GeV): Ecm_final= n/a"
/* << Ecm_final / 1_GeV*/
<<
endl
<<
"Elab_final="
<<
Elab_final
/
1
_GeV
<<
", Plab_final="
<<
(
Plab_final
/
1
_GeV
).
GetComponents
()
<<
", N_wounded,targ="
<<
QGSJetIIFragmentsStackData
::
GetWoundedNucleonsTarget
()
<<
", N_wounded,proj="
<<
QGSJetIIFragmentsStackData
::
GetWoundedNucleonsProjectile
()
<<
", N_fragm,proj="
<<
qfs
.
GetSize
()
<<
endl
;
}
return
process
::
EProcessReturn
::
eOk
;
}
}
// namespace corsika::process::qgsjetII
This diff is collapsed.
Click to expand it.
Processes/QGSJetII/Interaction.h
+
4
−
4
View file @
bc34c375
...
...
@@ -39,10 +39,10 @@ namespace corsika::process::qgsjetII {
corsika
::
particles
::
IsNucleus
(
TargetId
);
}
corsika
::
units
::
si
::
CrossSectionType
GetCrossSection
(
const
corsika
::
particles
::
Code
,
const
corsika
::
particles
::
Code
,
const
corsika
::
units
::
si
::
HEPEnergyType
,
const
unsigned
int
Abeam
=
0
,
const
unsigned
int
Atarget
=
0
)
const
;
corsika
::
units
::
si
::
CrossSectionType
GetCrossSection
(
const
corsika
::
particles
::
Code
,
const
corsika
::
particles
::
Code
,
const
corsika
::
units
::
si
::
HEPEnergyType
,
const
unsigned
int
Abeam
=
0
,
const
unsigned
int
Atarget
=
0
)
const
;
template
<
typename
TParticle
>
corsika
::
units
::
si
::
GrammageType
GetInteractionLength
(
TParticle
const
&
)
const
;
...
...
This diff is collapsed.
Click to expand it.
Preview
0%
Loading
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Save comment
Cancel
Please
register
or
sign in
to comment