diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index f6d04cfaf745817614ec9ae200b58d9591781e2b..101aa15ee56dd8711044fe1fcf42e812c96f1454 100644 --- a/Documentation/Examples/cascade_example.cc +++ b/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; diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h index 93516c141ad13888e47b57d72d247c897c2ebdbf..19f270eb2c80189d76a43b5f01cc85311eaad133 100644 --- a/Processes/Sibyll/Decay.h +++ b/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()));