use central boost and clean up

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;