IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 82e06a9b authored by Felix Riehn's avatar Felix Riehn
Browse files

adjust argument names in pythia

parent 4c838dec
No related branches found
No related tags found
No related merge requests found
......@@ -76,37 +76,37 @@ namespace corsika::process::pythia {
}
void Interaction::SetParticleListStable(
std::vector<particles::Code> const& particleList) {
for (auto p : particleList) Interaction::SetStable(p);
std::vector<particles::Code> const& vParticleList) {
for (auto p : vParticleList) Interaction::SetStable(p);
}
void Interaction::SetUnstable(const particles::Code pCode) {
cout << "Pythia::Interaction: setting " << pCode << " unstable.." << endl;
fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(pCode)), true);
void Interaction::SetUnstable(const particles::Code vCode) {
cout << "Pythia::Interaction: setting " << vCode << " unstable.." << endl;
fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(vCode)), true);
}
void Interaction::SetStable(const particles::Code pCode) {
cout << "Pythia::Interaction: setting " << pCode << " stable.." << endl;
fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(pCode)), false);
void Interaction::SetStable(const particles::Code vCode) {
cout << "Pythia::Interaction: setting " << vCode << " stable.." << endl;
fPythia.particleData.mayDecay(static_cast<int>(particles::GetPDG(vCode)), false);
}
void Interaction::ConfigureLabFrameCollision(
const particles::Code BeamId, const particles::Code TargetId,
const units::si::HEPEnergyType BeamEnergy) {
const particles::Code vBeamId, const particles::Code vTargetId,
const units::si::HEPEnergyType vBeamEnergy) {
using namespace units::si;
// Pythia configuration of the current event
// very clumsy. I am sure this can be done better..
// set beam
// beam id for pythia
auto const pdgBeam = static_cast<int>(particles::GetPDG(BeamId));
auto const pdgBeam = static_cast<int>(particles::GetPDG(vBeamId));
std::stringstream stBeam;
stBeam << "Beams:idA = " << pdgBeam;
fPythia.readString(stBeam.str());
// set target
auto pdgTarget = static_cast<int>(particles::GetPDG(TargetId));
auto pdgTarget = static_cast<int>(particles::GetPDG(vTargetId));
// replace hydrogen with proton, otherwise pythia goes into heavy ion mode!
if (TargetId == particles::Code::Hydrogen)
if (vTargetId == particles::Code::Hydrogen)
pdgTarget = static_cast<int>(particles::GetPDG(particles::Code::Proton));
std::stringstream stTarget;
stTarget << "Beams:idB = " << pdgTarget;
......@@ -114,7 +114,7 @@ namespace corsika::process::pythia {
// set frame to lab. frame
fPythia.readString("Beams:frameType = 2");
// set beam energy
const double Elab = BeamEnergy / 1_GeV;
const double Elab = vBeamEnergy / 1_GeV;
std::stringstream stEnergy;
stEnergy << "Beams:eA = " << Elab;
fPythia.readString(stEnergy.str());
......@@ -124,28 +124,28 @@ namespace corsika::process::pythia {
fPythia.init();
}
bool Interaction::CanInteract(const corsika::particles::Code pCode) {
return pCode == corsika::particles::Code::Proton ||
pCode == corsika::particles::Code::Neutron ||
pCode == corsika::particles::Code::AntiProton ||
pCode == corsika::particles::Code::AntiNeutron ||
pCode == corsika::particles::Code::PiMinus ||
pCode == corsika::particles::Code::PiPlus;
bool Interaction::CanInteract(const corsika::particles::Code vCode) {
return vCode == corsika::particles::Code::Proton ||
vCode == corsika::particles::Code::Neutron ||
vCode == corsika::particles::Code::AntiProton ||
vCode == corsika::particles::Code::AntiNeutron ||
vCode == corsika::particles::Code::PiMinus ||
vCode == corsika::particles::Code::PiPlus;
}
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
Interaction::GetCrossSection(const particles::Code BeamId,
const particles::Code TargetId,
const units::si::HEPEnergyType CoMenergy) {
Interaction::GetCrossSection(const particles::Code vBeamId,
const particles::Code vTargetId,
const units::si::HEPEnergyType vCoMenergy) {
using namespace units::si;
// interaction possible in pythia?
if (TargetId == particles::Code::Proton || TargetId == particles::Code::Hydrogen) {
if (CanInteract(BeamId) && ValidCoMEnergy(CoMenergy)) {
if (CanInteract(vBeamId) && ValidCoMEnergy(vCoMenergy)) {
// input particle PDG
auto const pdgCodeBeam = static_cast<int>(particles::GetPDG(BeamId));
auto const pdgCodeTarget = static_cast<int>(particles::GetPDG(TargetId));
const double ecm = CoMenergy / 1_GeV;
auto const pdgCodeBeam = static_cast<int>(particles::GetPDG(vBeamId));
auto const pdgCodeTarget = static_cast<int>(particles::GetPDG(vTargetId));
const double ecm = vCoMenergy / 1_GeV;
// calculate cross section
fSigma.calc(pdgCodeBeam, pdgCodeTarget, ecm);
......@@ -168,7 +168,7 @@ namespace corsika::process::pythia {
}
template <>
units::si::GrammageType Interaction::GetInteractionLength(Particle& p, Track&) {
units::si::GrammageType Interaction::GetInteractionLength(Particle& vP, Track&) {
using namespace units;
using namespace units::si;
......@@ -178,7 +178,7 @@ namespace corsika::process::pythia {
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
const particles::Code corsikaBeamId = p.GetPID();
const particles::Code corsikaBeamId = vP.GetPID();
// beam particles for pythia : 1, 2, 3 for p, pi, k
// read from cross section code table
......@@ -188,9 +188,9 @@ namespace corsika::process::pythia {
process::pythia::MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
// total momentum and energy
HEPEnergyType Elab = p.GetEnergy() + constants::nucleonMass;
HEPEnergyType Elab = vP.GetEnergy() + constants::nucleonMass;
process::pythia::MomentumVector pTotLab(rootCS, {0_GeV, 0_GeV, 0_GeV});
pTotLab += p.GetMomentum();
pTotLab += vP.GetMomentum();
pTotLab += pTarget;
auto const pTotLabNorm = pTotLab.norm();
// calculate cm. energy
......@@ -198,9 +198,9 @@ namespace corsika::process::pythia {
(Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
cout << "Interaction: LambdaInt: \n"
<< " input energy: " << p.GetEnergy() / 1_GeV << endl
<< " input energy: " << vP.GetEnergy() / 1_GeV << endl
<< " beam can interact:" << kInteraction << endl
<< " beam pid:" << p.GetPID() << endl;
<< " beam pid:" << vP.GetPID() << endl;
// TODO: move limits into variables
if (kInteraction && Elab >= 8.5_GeV && ValidCoMEnergy(ECoM)) {
......@@ -211,7 +211,7 @@ namespace corsika::process::pythia {
ideally as full particle object so that the four momenta
and the boosts can be defined..
*/
const auto* currentNode = p.GetNode();
const auto* currentNode = vP.GetNode();
const auto mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// determine average interaction length
......
......@@ -37,9 +37,9 @@ namespace corsika::process::pythia {
void SetStable(const corsika::particles::Code);
bool WasInitialized() { return fInitialized; }
bool ValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) {
bool ValidCoMEnergy(corsika::units::si::HEPEnergyType vEcm) {
using namespace corsika::units::si;
return (10_GeV < ecm) && (ecm < 1_PeV);
return (10_GeV < vEcm) && (vEcm < 1_PeV);
}
bool CanInteract(const corsika::particles::Code);
......@@ -47,9 +47,9 @@ namespace corsika::process::pythia {
const corsika::particles::Code,
const corsika::units::si::HEPEnergyType);
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
GetCrossSection(const corsika::particles::Code BeamId,
const corsika::particles::Code TargetId,
const corsika::units::si::HEPEnergyType CoMenergy);
GetCrossSection(const corsika::particles::Code,
const corsika::particles::Code,
const corsika::units::si::HEPEnergyType);
template <typename TParticle, typename TTrack>
corsika::units::si::GrammageType GetInteractionLength(TParticle&, TTrack&);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment