From 231c3d8cd333aaa7f2c009f4fb983cc971dc1b88 Mon Sep 17 00:00:00 2001 From: ralfulrich <ralf.ulrich@kit.edu> Date: Thu, 27 Dec 2018 09:49:33 +0100 Subject: [PATCH] temporarily fixed creation of sibyylCS in Interaction.h, ->100TeV shower at 45deg --- Documentation/Examples/cascade_example.cc | 6 +++--- Processes/Sibyll/Decay.h | 25 ++++++++++++++--------- Processes/Sibyll/Interaction.h | 3 ++- Processes/Sibyll/SibStack.h | 2 +- Processes/Sibyll/sibyll2.3c.cc | 2 +- Processes/Sibyll/sibyll2.3c.h | 4 ++-- 6 files changed, 24 insertions(+), 18 deletions(-) diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc index c3e618da1..451a3800d 100644 --- a/Documentation/Examples/cascade_example.cc +++ b/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(); diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h index 31609e426..74bc02418 100644 --- a/Processes/Sibyll/Decay.h +++ b/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()); diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h index 5dede2b41..e9edec807 100644 --- a/Processes/Sibyll/Interaction.h +++ b/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, diff --git a/Processes/Sibyll/SibStack.h b/Processes/Sibyll/SibStack.h index 8904ddaf2..bc3dc4c0f 100644 --- a/Processes/Sibyll/SibStack.h +++ b/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; } diff --git a/Processes/Sibyll/sibyll2.3c.cc b/Processes/Sibyll/sibyll2.3c.cc index 775c97e66..f614b596b 100644 --- a/Processes/Sibyll/sibyll2.3c.cc +++ b/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 = diff --git a/Processes/Sibyll/sibyll2.3c.h b/Processes/Sibyll/sibyll2.3c.h index 00a1cfa95..104b55eb0 100644 --- a/Processes/Sibyll/sibyll2.3c.h +++ b/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&); } -- GitLab