From 15227373e8de74c94e72697eb19017856dac6f09 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Tue, 25 Dec 2018 13:39:08 +0100
Subject: [PATCH] TrackWriter output to ofstream
---
Documentation/Examples/cascade_example.cc | 5 +-
Framework/Cascade/Cascade.h | 20 +++----
Processes/Sibyll/Interaction.h | 65 +++++++++++-----------
Processes/TrackWriter/TrackWriter.h | 4 +-
Processes/TrackingLine/TrackingLine.h | 2 +-
Processes/TrackingLine/testTrackingLine.cc | 2 +
6 files changed, 51 insertions(+), 47 deletions(-)
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 0e6ce843b..323617017 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -7,7 +7,7 @@
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
-
+
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/stack_inspector/StackInspector.h>
@@ -191,7 +191,6 @@ public:
// The example main program for a particle cascade
//
int main() {
-
// initialize random number sequence(s)
corsika::random::RNGManager::GetInstance().RegisterRandomStream("cascade");
@@ -229,7 +228,7 @@ int main() {
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
- const hep::EnergyType E0 = 100_TeV;
+ const hep::EnergyType E0 = 10_TeV;
{
auto particle = stack.NewParticle();
particle.SetPID(Code::Proton);
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index fc0dfffb9..6aed37d64 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -57,11 +57,7 @@ namespace corsika::cascade {
}
void Step(Particle& particle) {
-
using namespace corsika::units::si;
- using std::cout;
- using std::endl;
- using std::log;
// determine geometric tracking
corsika::setup::Trajectory step = fTracking.GetTrack(particle);
@@ -78,11 +74,15 @@ namespace corsika::cascade {
<< ", next_interact=" << next_interact << std::endl;
// convert next_step from grammage to length
+ auto const* currentNode =
+ fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
+
+ if (currentNode == &*fEnvironment.GetUniverse()) {
+ throw std::runtime_error("particle entered void universe");
+ }
+
LengthType const distance_interact =
- fEnvironment.GetUniverse()
- ->GetContainingNode(particle.GetPosition())
- ->GetModelProperties()
- .ArclengthFromGrammage(step, next_interact);
+ currentNode->GetModelProperties().ArclengthFromGrammage(step, next_interact);
// determine the maximum geometric step length
LengthType const distance_max = fProcessSequence.MaxStepLength(particle, step);
@@ -113,7 +113,7 @@ namespace corsika::cascade {
// std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
particle.SetPosition(step.PositionFromArclength(min_distance));
// .... also update time, momentum, direction, ...
-
+
step.LimitEndTo(min_distance);
// apply all continuous processes on particle + track
@@ -127,7 +127,7 @@ namespace corsika::cascade {
return;
}
- std::cout << "sth. happening before geometric limit ?"
+ std::cout << "sth. happening before geometric limit ? "
<< ((min_distance < distance_max) ? "yes" : "no") << std::endl;
if (min_distance < distance_max) { // interaction to happen within geometric limit
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 23af2bf2a..fe8cfecef 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -64,7 +64,8 @@ namespace corsika::process::sibyll {
// FOR NOW: assume target is oxygen
const int kTarget = corsika::particles::Oxygen::GetNucleusA();
- const hep::MassType nucleon_mass = 0.5*(corsika::particles::Proton::GetMass()+corsika::particles::Neutron::GetMass());
+ const hep::MassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
+ corsika::particles::Neutron::GetMass());
hep::EnergyType Etot = p.GetEnergy() + nucleon_mass;
MomentumVector Ptot(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
// FOR NOW: assume target is at rest
@@ -128,22 +129,23 @@ namespace corsika::process::sibyll {
<< process::sibyll::CanInteract(p.GetPID()) << endl;
if (process::sibyll::CanInteract(p.GetPID())) {
const CoordinateSystem& rootCS =
- RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
Point pOrig = p.GetPosition();
TimeType tOrig = p.GetTime();
- // sibyll CS has z along particle momentum
- // FOR NOW: hard coded z-axis for corsika frame
- QuantityVector<length_d> const zAxis{0_m, 0_m, 1_m};
- QuantityVector<length_d> const xAxis{1_m, 0_m, 0_m};
- auto pt = []( MomentumVector p ){
- return sqrt(p.GetComponents()[0] * p.GetComponents()[0] + p.GetComponents()[1] * p.GetComponents()[1]);
- };
- double theta = acos( p.GetMomentum().GetComponents()[2] / p.GetMomentum().norm());
- cout << "ProcessSibyll: polar angle between sibyllCS and rootCS: " << theta << endl;
- double phi = asin( p.GetMomentum().GetComponents()[0]/pt(p.GetMomentum() ) );
- CoordinateSystem sibyllCS = rootCS.rotate(xAxis, theta);
-
+ // sibyll CS has z along particle momentum
+ // FOR NOW: hard coded z-axis for corsika frame
+ QuantityVector<length_d> const zAxis{0_m, 0_m, 1_m};
+ QuantityVector<length_d> const xAxis{1_m, 0_m, 0_m};
+ auto pt = [](MomentumVector p) {
+ return sqrt(p.GetComponents()[0] * p.GetComponents()[0] +
+ p.GetComponents()[1] * p.GetComponents()[1]);
+ };
+ double theta = acos(p.GetMomentum().GetComponents()[2] / p.GetMomentum().norm());
+ cout << "ProcessSibyll: polar angle between sibyllCS and rootCS: " << theta
+ << endl;
+ CoordinateSystem sibyllCS = rootCS.rotate(xAxis, theta);
+
/*
the target should be defined by the Environment,
ideally as full particle object so that the four momenta
@@ -154,17 +156,18 @@ namespace corsika::process::sibyll {
*/
// FOR NOW: set target to oxygen
const int kTarget = corsika::particles::Oxygen::
- GetNucleusA(); // env.GetTargetParticle().GetPID();
+ GetNucleusA(); // env.GetTargetParticle().GetPID();
// FOR NOW: target is always at rest
- const hep::MassType nucleon_mass = 0.5*(corsika::particles::Proton::GetMass()+corsika::particles::Neutron::GetMass());
+ const hep::MassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
+ corsika::particles::Neutron::GetMass());
const EnergyType Etarget = 0_GeV + nucleon_mass;
const auto pTarget = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
cout << "target momentum (GeV/c): " << pTarget.GetComponents() / 1_GeV << endl;
cout << "beam momentum (GeV/c): " << p.GetMomentum().GetComponents() / 1_GeV
<< endl;
- cout << "beam momentum in sibyll frame (GeV/c): " << p.GetMomentum().GetComponents(sibyllCS) / 1_GeV
- << endl;
+ cout << "beam momentum in sibyll frame (GeV/c): "
+ << p.GetMomentum().GetComponents(sibyllCS) / 1_GeV << endl;
cout << "position of interaction: " << pOrig.GetCoordinates() << endl;
cout << "time: " << tOrig << endl;
@@ -224,31 +227,31 @@ namespace corsika::process::sibyll {
// add particles from sibyll to stack
// link to sibyll stack
- // here we need to pass projectile momentum and energy to define the local sibyll frame
- // and the boosts to the lab. frame
+ // here we need to pass projectile momentum and energy to define the local
+ // sibyll frame and the boosts to the lab. frame
SibStack ss;
// momentum array in Sibyll
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
EnergyType E_final = 0_GeV, Ecm_final = 0_GeV;
for (auto& psib : ss) {
-
+
// skip particles that have decayed in Sibyll
- if( psib.HasDecayed()) continue;
+ if (psib.HasDecayed()) continue;
// transform energy to lab. frame, primitve
// compute beta_vec * p_vec
// arbitrary Lorentz transformation based on sibyll routines
const auto gammaBetaComponents = gambet.GetComponents();
- // FOR NOW: fill vector in sibCS and then rotate into rootCS
- // can be done in SibStack by passing sibCS
- // get momentum vector in sibyllCS
+ // FOR NOW: fill vector in sibCS and then rotate into rootCS
+ // can be done in SibStack by passing sibCS
+ // get momentum vector in sibyllCS
const auto pSibyllComponentsSibCS = psib.GetMomentum().GetComponents();
- // temporary vector in sibyllCS
- auto SibVector = MomentumVector( sibyllCS, pSibyllComponentsSibCS);
- // rotatate to rootCS
- const auto pSibyllComponents = SibVector.GetComponents(rootCS);
- // boost to lab. frame
+ // temporary vector in sibyllCS
+ auto SibVector = MomentumVector(sibyllCS, pSibyllComponentsSibCS);
+ // rotatate to rootCS
+ const auto pSibyllComponents = SibVector.GetComponents(rootCS);
+ // boost to lab. frame
hep::EnergyType en_lab = 0. * 1_GeV;
hep::MomentumType p_lab_components[3];
en_lab = psib.GetEnergy() * gamma;
diff --git a/Processes/TrackWriter/TrackWriter.h b/Processes/TrackWriter/TrackWriter.h
index b86ae5ed3..50e17e6af 100644
--- a/Processes/TrackWriter/TrackWriter.h
+++ b/Processes/TrackWriter/TrackWriter.h
@@ -36,7 +36,7 @@ namespace corsika::process::TrackWriter {
auto const delta = t.GetPosition(1).GetCoordinates() - start;
auto const& name = corsika::particles::GetName(p.GetPID());
- std::cerr << name << " " << start[0] / 1_m << ' ' << start[1] / 1_m << ' '
+ fFile << name << " " << 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';
@@ -50,7 +50,7 @@ namespace corsika::process::TrackWriter {
private:
std::string const fFilename;
- std::ofstream fFile;
+ mutable std::ofstream fFile;
};
} // namespace corsika::process::TrackWriter
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 6fdd4f662..6b554feba 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -84,7 +84,7 @@ namespace corsika::process {
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()
<< std::endl;
std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
index 737b536f2..285f275c0 100644
--- a/Processes/TrackingLine/testTrackingLine.cc
+++ b/Processes/TrackingLine/testTrackingLine.cc
@@ -9,6 +9,7 @@
*/
#include <corsika/environment/Environment.h>
+#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/tracking_line/TrackingLine.h>
@@ -47,6 +48,7 @@ struct DummyParticle {
auto GetEnergy() const { return fEnergy; }
auto GetMomentum() const { return fMomentum; }
auto GetPosition() const { return fPosition; }
+ auto GetPID() const { return corsika::particles::Code::Unknown; }
};
struct DummyStack {
--
GitLab