diff --git a/Processes/Pythia/Decay.cc b/Processes/Pythia/Decay.cc index 27890d28bfe1ce42576fc81d08550ea80daf70de..a17678d18c9932a1f6417f792ab5cb0e40946f6a 100644 --- a/Processes/Pythia/Decay.cc +++ b/Processes/Pythia/Decay.cc @@ -57,7 +57,7 @@ namespace corsika::process::pythia { fPythia.readString("Next:numberShowEvent = 0"); fPythia.readString("Print:quiet = on"); - fPythia.readString("Check:particleData = 1"); + fPythia.readString("Check:particleData = 0"); /* switching off event check in pythia is needed to allow decays that are off-shell @@ -65,7 +65,7 @@ namespace corsika::process::pythia { the consistency of particle masses between event generators is an unsolved issues */ cout << "Pythia::Init: switching off event checking in pythia.." << endl; - fPythia.readString("Check:event = 0"); + fPythia.readString("Check:event = 1"); fPythia.readString("ProcessLevel:all = off"); fPythia.readString("ProcessLevel:resonanceDecays = off"); diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc index 916ea2ec0829f6de41b4055428bac17b78f0fb62..e7d610dc59831d026503583fb49402f1e0eaffb2 100644 --- a/Processes/Sibyll/Interaction.cc +++ b/Processes/Sibyll/Interaction.cc @@ -316,20 +316,50 @@ namespace corsika::process::sibyll { HEPEnergyType const eCoM = psib.GetEnergy(); auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM)); + auto const pid = process::sibyll::ConvertFromSibyll(psib.GetPID()); + // check if on-shell in corsika + auto const m_kinetic = Plab.GetNorm(); + auto const m_corsika = particles::GetMass(pid); + auto const e_corsika = Plab.GetTimeLikeComponent(); + auto const m_sibyll = corsika::process::sibyll::GetSibyllMass(pid); + auto const m_err = abs(m_kinetic - m_corsika) / m_corsika; + if (m_err > 1.e-5) { + const HEPEnergyType e_shift_corsika = sqrt( + Plab.GetSpaceLikeComponents().GetSquaredNorm() + m_corsika * m_corsika); + auto const e_shift_relative = (e_shift_corsika / e_corsika - 1) * 100; + // warn about percent level shifts in particle energy + if (abs(e_shift_relative) > 1) { + std::cout << "Sibyll::Interaction: shifted particle energy by " + << e_shift_relative << " %" << std::endl; + + std::cout << "shift particle mass for " << pid << std::endl + << "corsika mass (GeV): " << m_corsika / 1_GeV << std::endl + << "kinetic mass (GeV): " << m_kinetic / 1_GeV << std::endl + << "sibyll mass (GeV): " << m_sibyll / 1_GeV << std::endl + << "(m_kin-m_cor)/en: " << m_err << std::endl; + } + Plab.GetTimeLikeComponent() = e_shift_corsika; + } + // add to corsika stack auto pnew = vP.AddSecondary( tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector, geometry::Point, units::si::TimeType>{ - process::sibyll::ConvertFromSibyll(psib.GetPID()), - Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig, + pid, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig, tOrig}); Plab_final += pnew.GetMomentum(); Elab_final += pnew.GetEnergy(); Ecm_final += psib.GetEnergy(); } - cout << "conservation (all GeV): Ecm_final=" << Ecm_final / 1_GeV << endl - << "Elab_final=" << Elab_final / 1_GeV + cout << "conservation (all GeV):" << endl + << "Ecm_initial=" << Ecm / 1_GeV << " Ecm_final=" << Ecm_final / 1_GeV + << endl + << "Elab_initial=" << eProjectileLab / 1_GeV + << " Elab_final=" << Elab_final / 1_GeV + << " diff (%)=" << (Elab_final / eProjectileLab / get_nwounded() - 1) * 100 + << endl + << "Plab_initial=" << (pProjectileLab / 1_GeV).GetComponents() << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl; } }