temporarily fixed creation of sibyylCS in Interaction.h, ->100TeV shower at 45deg

Documentation/Examples/cascade_example.cc

@@ -80,7 +80,7 @@ public:
       case Code::Electron:
         is_em = true;
         break;
-    case Code::Positron:
+      case Code::Positron:
         is_em = true;
         break;
       case Code::Gamma:
@@ -237,8 +237,8 @@ int main() {
   // setup particle stack, and add primary particle
   setup::Stack stack;
   stack.Clear();
-  const hep::EnergyType E0 = 10_TeV;
-  double theta = 0.;
+  const hep::EnergyType E0 = 100_TeV;
+  double theta = 45.;
   double phi = 0.;
   {
     auto particle = stack.NewParticle();

Processes/Sibyll/Decay.h

@@ -145,15 +145,19 @@ namespace corsika::process {
             corsika::particles::GetLifetime(p.GetPID());
         auto const lifetime = gamma * t0;
 
-	const auto mkin =  (E * E - p.GetMomentum().squaredNorm());//delta_mass(p.GetMomentum(), E, m);
+        const auto mkin =
+            (E * E - p.GetMomentum().squaredNorm()); // delta_mass(p.GetMomentum(), E, m);
         cout << "Decay: code: " << p.GetPID() << endl;
         cout << "Decay: MinStep: t0: " << t0 << endl;
         cout << "Decay: MinStep: energy: " << E / 1_GeV << " GeV" << endl;
-	cout << "Decay: momentum: " <<  p.GetMomentum().GetComponents() / 1_GeV << " GeV" << endl;
-	cout << "Decay: momentum: shell mass-kin. inv. mass " <<  mkin / 1_GeV / 1_GeV << " " << m / 1_GeV*m / 1_GeV << endl;
-	//cout << "Decay: sib mass: " << s_mass1_.am2[ process::sibyll::ConvertToSibyllRaw(p.GetPID()) ] << endl;
-	auto sib_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
-	cout << "Decay: sib mass: " << get_sibyll_mass2( sib_id )  << endl;
+        cout << "Decay: momentum: " << p.GetMomentum().GetComponents() / 1_GeV << " GeV"
+             << endl;
+        cout << "Decay: momentum: shell mass-kin. inv. mass " << mkin / 1_GeV / 1_GeV
+             << " " << m / 1_GeV * m / 1_GeV << endl;
+        // cout << "Decay: sib mass: " << s_mass1_.am2[
+        // process::sibyll::ConvertToSibyllRaw(p.GetPID()) ] << endl;
+        auto sib_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+        cout << "Decay: sib mass: " << get_sibyll_mass2(sib_id) << endl;
         cout << "Decay: MinStep: gamma: " << gamma << endl;
         cout << "Decay: MinStep: tau: " << lifetime << endl;
 
@@ -176,13 +180,14 @@ namespace corsika::process {
         pin.SetPID(process::sibyll::ConvertToSibyllRaw(pCode));
         pin.SetEnergy(p.GetEnergy());
         pin.SetMomentum(p.GetMomentum());
-	// setting particle mass with Corsika values, may be inconsistent with sibyll internal values
+        // setting particle mass with Corsika values, may be inconsistent with sibyll
+        // internal values
 #warning setting particle mass with Corsika values, may be inconsistent with sibyll internal values
-	pin.SetMass( corsika::particles::GetMass( pCode ) );
+        pin.SetMass(corsika::particles::GetMass(pCode));
         // remember position
         Point decayPoint = p.GetPosition();
         TimeType t0 = p.GetTime();
-	// remove original particle from corsika stack
+        // remove original particle from corsika stack
         p.Delete();
         // set all particles/hadrons unstable
         // setHadronsUnstable();
@@ -199,7 +204,7 @@ namespace corsika::process {
         // copy particles from sibyll stack to corsika
         for (auto& psib : ss) {
           // FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
-          if ( psib.HasDecayed() ) continue;
+          if (psib.HasDecayed()) continue;
           // add to corsika stack
           auto pnew = s.NewParticle();
           pnew.SetEnergy(psib.GetEnergy());

Processes/Sibyll/Interaction.h

@@ -171,7 +171,8 @@ namespace corsika::process::sibyll {
         cout << "ProcessSibyll: azimuth angle between sibyllCS and rootCS: "
              << phi / M_PI * 180. << endl;
         // double phi = asin( p.GetMomentum().GetComponents()[0]/pt(p.GetMomentum() ) );
-        CoordinateSystem sibyllCS = rootCS.rotate(zAxis, phi).rotate(yAxis, theta);
+        const CoordinateSystem tempCS = rootCS.rotate(zAxis, phi);
+        const CoordinateSystem sibyllCS = tempCS.rotate(yAxis, theta);
 
         /*
            the target should be defined by the Environment,

Processes/Sibyll/SibStack.h

@@ -107,7 +107,7 @@ namespace corsika::process::sibyll {
     corsika::units::hep::EnergyType GetMass() const {
       return GetStackData().GetMass(GetIndex());
     }
-   
+
     bool HasDecayed() const {
       return abs(GetStackData().GetId(GetIndex())) > 100 ? true : false;
     }

Processes/Sibyll/sibyll2.3c.cc

@@ -15,7 +15,7 @@
 #include <random>
 
 int get_nwounded() { return s_chist_.nwd; }
-double get_sibyll_mass2( int& id ) { return  s_mass1_.am2[ id ]; }
+double get_sibyll_mass2(int& id) { return s_mass1_.am2[id]; }
 
 double s_rndm_(int&) {
   static corsika::random::RNG& rng =

Processes/Sibyll/sibyll2.3c.h

@@ -100,9 +100,9 @@ void sib_sigma_hp_(const int&, const double&, double&, double&, double&, double*
 
 double s_rndm_(int&);
 
-int get_nwounded();  
+int get_nwounded();
 double get_sibyll_mass2(int&);
-  
+
 // phojet random generator setup
 void pho_rndin_(int&, int&, int&, int&);
 }