TrackWriter output to ofstream

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);

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

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;

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

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;

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 {