repaired interface to sibyll decay routine to properly pass particle masses,...

repaired interface to sibyll decay routine to properly pass particle masses, added positrons as EM particle to EMCut

Documentation/Examples/cascade_example.cc

@@ -75,6 +75,9 @@ public:
       case Code::Electron:
         is_em = true;
         break;
+    case Code::Positron:
+        is_em = true;
+        break;
       case Code::Gamma:
         is_em = true;
         break;

Processes/Sibyll/Decay.h

@@ -132,9 +132,15 @@ 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);
         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: MinStep: gamma: " << gamma << endl;
         cout << "Decay: MinStep: tau: " << lifetime << endl;
 
@@ -154,10 +160,13 @@ 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
+#warning setting particle mass with Corsika values, may be inconsistent with sibyll internal values
+	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();
@@ -170,15 +179,12 @@ namespace corsika::process {
         // print output
         int print_unit = 6;
         sib_list_(print_unit);
+
         // copy particles from sibyll stack to corsika
-        int i = -1;
         for (auto& psib : ss) {
-          ++i;
           // FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
-          if (abs(s_plist_.llist[i]) > 100) continue;
+          if ( psib.HasDecayed() ) continue;
           // add to corsika stack
-          // cout << "decay product: " << process::sibyll::ConvertFromSibyll(
-          // psib.GetPID() ) << endl;
           auto pnew = s.NewParticle();
           pnew.SetEnergy(psib.GetEnergy());
           pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));