From 5074f107b135c3cbcf37d879936dc1c6e69bd93a Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Sat, 12 Jan 2019 18:53:33 +0000
Subject: [PATCH] use central boost and clean up
---
Processes/Sibyll/Interaction.h | 188 +++++++++++----------------------
1 file changed, 60 insertions(+), 128 deletions(-)
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index b8709216..e646858e 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -178,66 +178,55 @@ namespace corsika::process::sibyll {
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ // position and time of interaction, not used in Sibyll
Point pOrig = p.GetPosition();
TimeType tOrig = p.GetTime();
- // sibyll CS has z along particle momentum
- // FOR NOW: hard coded z-axis for corsika frame
- QuantityVector<length_d> const zAxis{0_m, 0_m, 1_m};
- QuantityVector<length_d> const yAxis{0_m, 1_m, 0_m};
- auto rotation_angles = [](MomentumVector const& pin) {
- const auto p = pin.GetComponents();
- const auto th = acos(p[2] / p.norm());
- const auto ph = atan2(
- p[1] / 1_GeV, p[0] / 1_GeV); // acos( p[0] / sqrt(p[0]*p[0]+p[1]*p[1] ) );
- return std::make_tuple(th, ph);
- };
- // auto pt = []( MomentumVector &p ){
- // return sqrt(p.GetComponents()[0] * p.GetComponents()[0] +
- // p.GetComponents()[1] * p.GetComponents()[1]);
- // };
- // double theta = acos( p.GetMomentum().GetComponents()[2] /
- // p.GetMomentum().norm());
- auto const [theta, phi] = rotation_angles(p.GetMomentum());
- cout << "ProcessSibyll: zenith angle between sibyllCS and rootCS: "
- << theta / M_PI * 180. << endl;
- cout << "ProcessSibyll: azimuth angle between sibyllCS and rootCS: "
- << phi / M_PI * 180. << endl;
- // double phi = asin( p.GetMomentum().GetComponents()[0]/pt(p.GetMomentum() ) );
- const CoordinateSystem tempCS = rootCS.rotate(zAxis, phi);
- const CoordinateSystem sibyllCS = tempCS.rotate(yAxis, theta);
-
-
- // FOR NOW: target is always at rest
+ // define target
+ // for Sibyll is always a single nucleon
auto constexpr nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
corsika::particles::Neutron::GetMass());
- const auto Etarget = 0_GeV + nucleon_mass;
- const auto pTarget = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
- cout << "target momentum (GeV/c): " << pTarget.GetComponents() / 1_GeV << endl;
- cout << "beam momentum (GeV/c): " << p.GetMomentum().GetComponents() / 1_GeV
- << endl;
- cout << "beam momentum in sibyll frame (GeV/c): "
- << p.GetMomentum().GetComponents(sibyllCS) / 1_GeV << endl;
-
- cout << "position of interaction: " << pOrig.GetCoordinates() << endl;
- cout << "time: " << tOrig << endl;
-
- // get energy of particle from stack
- /*
- stack is in GeV in lab. frame
- convert to GeV in cm. frame
- (assuming proton at rest as target AND
- assuming no pT, i.e. shower frame-z is aligned with hadron-int-frame-z)
- */
- // total energy: E_beam + E_target
- // in lab. frame: E_beam + m_target*c**2
- HEPEnergyType E = p.GetEnergy();
- HEPEnergyType Etot = E + Etarget;
- // total momentum
- MomentumVector Ptot = p.GetMomentum();
- // invariant mass, i.e. cm. energy
- HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
+ // FOR NOW: target is always at rest
+ const auto eTargetLab = 0_GeV + nucleon_mass;
+ const auto pTargetLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
+
+ // define projectile
+ auto const pProjectileLab = p.GetMomentum();
+ HEPEnergyType const eProjectileLab = p.GetEnergy();
+
+ // define target kinematics in lab frame
+ HEPMassType const targetMass = nucleon_mass;
+ // define boost to and from CoM frame
+ // CoM frame definition in Sibyll projectile: +z
+ COMBoost const boost(eProjectileLab, pProjectileLab, targetMass);
+ cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << endl
+ << "Interaction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV << endl;
+ cout << "Interaction: etarget lab: " << eTargetLab / 1_GeV << endl
+ << "Interaction: ptarget lab: " << pTargetLab.GetComponents() / 1_GeV << endl;
+
+ // just for show:
+ // boost projecticle
+ auto const [eProjectileCoM, pProjectileCoM] =
+ boost.toCoM(eProjectileLab, pProjectileLab);
+
+ // boost target
+ auto const [eTargetCoM, pTargetCoM] =
+ boost.toCoM(eTargetLab, pTargetLab);
+
+ cout << "Interaction: ebeam CoM: " << eProjectileCoM / 1_GeV << endl
+ << "Interaction: pbeam CoM: " << pProjectileCoM / 1_GeV << endl;
+ cout << "Interaction: etarget CoM: " << eTargetCoM / 1_GeV << endl
+ << "Interaction: ptarget CoM: " << pTargetCoM / 1_GeV << endl;
+
+ cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl;
+ cout << "Interaction: time: " << tOrig << endl;
+
+ HEPEnergyType Etot = eProjectileLab + eTargetLab;
+ MomentumVector Ptot = p.GetMomentum();
+ // invariant mass, i.e. cm. energy
+ HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
+
// sample target mass number
const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
@@ -270,45 +259,13 @@ namespace corsika::process::sibyll {
if(targetSibCode>18||targetSibCode<1)
throw std::runtime_error("Sibyll target outside range. Only nuclei with A<18 or protons are allowed.");
- /*
- get transformation between Stack-frame and SibStack-frame
- for EAS Stack-frame is lab. frame, could be different for CRMC-mode
- the transformation should be derived from the input momenta
- */
- const double gamma = Etot / Ecm;
- const auto gambet = Ptot / Ecm;
-
- std::cout << "Interaction: "
- << " DoDiscrete: gamma:" << gamma << endl;
- std::cout << "Interaction: "
- << " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
-
- auto const pProjectileLab = p.GetMomentum();
- //{rootCS, {0_GeV / c, 1_PeV / c, 0_GeV / c}};
- HEPEnergyType const eProjectileLab = p.GetEnergy();
- //energy(projectileMass, pProjectileLab);
-
- // define target kinematics in lab frame
- HEPMassType const targetMass = nucleon_mass;
- // define boost to com frame
- COMBoost const boost(eProjectileLab, pProjectileLab, targetMass);
-
- cout << "Interaction: new boost: ebeam lab: " << eProjectileLab / 1_GeV << endl
- << "Interaction: new boost: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV << endl;
-
- // boost projecticle
- auto const [eProjectileCoM, pProjectileCoM] =
- boost.toCoM(eProjectileLab, pProjectileLab);
-
- cout << "Interaction: new boost: ebeam com: " << eProjectileCoM / 1_GeV << endl
- << "Interaction: new boost: pbeam com: " << pProjectileCoM / 1_GeV << endl;
-
+ // beam id for sibyll
const int kBeam = process::sibyll::ConvertToSibyllRaw(corsikaBeamId);
std::cout << "Interaction: "
- << " DoInteraction: E(GeV):" << E / 1_GeV
+ << " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
<< " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
- if (E < 8.5_GeV || Ecm < 10_GeV) {
+ if ( eProjectileLab < 8.5_GeV || Ecm < 10_GeV) {
std::cout << "Interaction: "
<< " DoInteraction: should have dropped particle.." << std::endl;
// p.Delete(); delete later... different process
@@ -331,62 +288,37 @@ namespace corsika::process::sibyll {
// add particles from sibyll to stack
// link to sibyll stack
- // here we need to pass projectile momentum and energy to define the local
- // sibyll frame and the boosts to the lab. frame
SibStack ss;
- // momentum array in Sibyll
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
- HEPEnergyType E_final = 0_GeV, Ecm_final = 0_GeV;
+ HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
for (auto& psib : ss) {
// skip particles that have decayed in Sibyll
if (psib.HasDecayed()) continue;
- // transform energy to lab. frame, primitve
- // compute beta_vec * p_vec
- // arbitrary Lorentz transformation based on sibyll routines
- const auto gammaBetaComponents = gambet.GetComponents();
- // FOR NOW: fill vector in sibCS and then rotate into rootCS
- // can be done in SibStack by passing sibCS
- // get momentum vector in sibyllCS
- const auto pSibyllComponentsSibCS = psib.GetMomentum().GetComponents();
- // temporary vector in sibyllCS
- auto SibVector = MomentumVector(sibyllCS, pSibyllComponentsSibCS);
- // rotatate to rootCS
- const auto pSibyllComponents = SibVector.GetComponents(rootCS);
- // boost to lab. frame
- HEPEnergyType en_lab = 0. * 1_GeV;
- HEPMomentumType p_lab_components[3];
- en_lab = psib.GetEnergy() * gamma;
- HEPEnergyType pnorm = 0. * 1_GeV;
- for (int j = 0; j < 3; ++j)
- pnorm += (pSibyllComponents[j] * gammaBetaComponents[j]) / (gamma + 1.);
- pnorm += psib.GetEnergy();
-
- for (int j = 0; j < 3; ++j) {
- p_lab_components[j] =
- pSibyllComponents[j] - (-1) * pnorm * gammaBetaComponents[j];
- en_lab -= (-1) * pSibyllComponents[j] * gammaBetaComponents[j];
- }
-
+ // transform energy to lab. frame
+ auto const pCoM = psib.GetMomentum().GetComponents();
+ HEPEnergyType const eCoM = psib.GetEnergy();
+ auto const [ eLab, pLab ] = boost.fromCoM( eCoM, pCoM );
+
// add to corsika stack
auto pnew = s.NewParticle();
- pnew.SetEnergy(en_lab);
- pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
- corsika::geometry::QuantityVector<hepmomentum_d> p_lab_c{
- p_lab_components[0], p_lab_components[1], p_lab_components[2]};
- pnew.SetMomentum(MomentumVector(rootCS, p_lab_c));
+ pnew.SetPID(process::sibyll::ConvertFromSibyll(psib.GetPID()));
+ pnew.SetEnergy(eLab);
+ pnew.SetMomentum(pLab);
pnew.SetPosition(pOrig);
pnew.SetTime(tOrig);
+
Plab_final += pnew.GetMomentum();
- E_final += en_lab;
+ Elab_final += pnew.GetEnergy();
Ecm_final += psib.GetEnergy();
}
- std::cout << "conservation (all GeV): E_final=" << E_final / 1_GeV
- << ", Ecm_final=" << Ecm_final / 1_GeV
+ std::cout << "conservation (all GeV): Ecm_final=" << Ecm_final / 1_GeV << std::endl
+ << "Elab_final=" << Elab_final / 1_GeV
<< ", Plab_final=" << (Plab_final / 1_GeV).GetComponents()
<< std::endl;
+
}
}
return process::EProcessReturn::eOk;
--
GitLab