IAP GITLAB

Commit 2f4918e1 authored by Andre Schmidt's avatar Andre Schmidt

save before change

parent f2e86927
Pipeline #2035 failed with stages
in 0 seconds
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <catch2/catch.hpp>
#include <corsika/geometry/FourVector.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/COMBoost.h>
#include <iostream>
using namespace corsika::geometry;
using namespace corsika::utl;
using namespace corsika::units::si;
using corsika::units::constants::c;
using corsika::units::constants::cSquared;
double constexpr absMargin = 1e-6;
CoordinateSystem const& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// helper function for energy-momentum
// relativistic energy
auto const energy = [](HEPMassType m, Vector<hepmomentum_d> const& p) {
return sqrt(m * m + p.squaredNorm());
};
// helper function for mandelstam-s
auto const s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) {
return E * E - p.squaredNorm();
};
TEST_CASE("boosts") {
// define target kinematics in lab frame
HEPMassType const targetMass = 1_GeV;
Vector<hepmomentum_d> pTargetLab{rootCS, {0_eV, 0_eV, 0_eV}};
HEPEnergyType const eTargetLab = energy(targetMass, pTargetLab);
/*
General tests check the interface and basic operation
*/
SECTION("General tests") {
// define projectile kinematics in lab frame
HEPMassType const projectileMass = 1._GeV;
Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 1_PeV, 0_GeV}};
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
const FourVector PprojLab(eProjectileLab, pProjectileLab);
// define boost to com frame
COMBoost boost(PprojLab, targetMass);
// boost projecticle
auto const PprojCoM = boost.toCoM(PprojLab);
// boost target
auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
// sum of momenta in CoM, should be 0
auto const sumPCoM =
PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
CHECK(sumPCoM.norm() / 1_GeV == Approx(0).margin(absMargin));
// mandelstam-s should be invariant under transformation
CHECK(s(eProjectileLab + eTargetLab,
pProjectileLab.GetComponents() + pTargetLab.GetComponents()) /
1_GeV / 1_GeV ==
Approx(s(PprojCoM.GetTimeLikeComponent() + PtargCoM.GetTimeLikeComponent(),
PprojCoM.GetSpaceLikeComponents().GetComponents() +
PtargCoM.GetSpaceLikeComponents().GetComponents()) /
1_GeV / 1_GeV));
// boost back...
auto const PprojBack = boost.fromCoM(PprojCoM);
// ...should yield original values before the boosts
CHECK(PprojBack.GetTimeLikeComponent() / PprojLab.GetTimeLikeComponent() ==
Approx(1));
CHECK(
(PprojBack.GetSpaceLikeComponents() - PprojLab.GetSpaceLikeComponents()).norm() /
PprojLab.GetSpaceLikeComponents().norm() ==
Approx(0).margin(absMargin));
}
/*
special case: projectile along -z
*/
SECTION("Test boost along z-axis") {
// define projectile kinematics in lab frame
HEPMassType const projectileMass = 1_GeV;
Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_PeV, -1_PeV}};
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
const FourVector PprojLab(eProjectileLab, pProjectileLab);
// define boost to com frame
COMBoost boost(PprojLab, targetMass);
// boost projecticle
auto const PprojCoM = boost.toCoM(PprojLab);
// boost target
auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
// sum of momenta in CoM, should be 0
auto const sumPCoM =
PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
CHECK(sumPCoM.norm() / 1_GeV == Approx(0).margin(absMargin));
}
/*
special case: projectile with arbitrary direction
*/
SECTION("Test boost along tilted axis") {
const HEPMomentumType P0 = 1_PeV;
double theta = 33.;
double phi = -10.;
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
};
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
// define projectile kinematics in lab frame
HEPMassType const projectileMass = 1_GeV;
Vector<hepmomentum_d> pProjectileLab(rootCS, {px, py, pz});
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
const FourVector PprojLab(eProjectileLab, pProjectileLab);
// define boost to com frame
COMBoost boost(PprojLab, targetMass);
// boost projecticle
auto const PprojCoM = boost.toCoM(PprojLab);
// boost target
auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
// sum of momenta in CoM, should be 0
auto const sumPCoM =
PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
CHECK(sumPCoM.norm() / 1_GeV == Approx(0).margin(absMargin));
}
/*
test the ultra-high energy behaviour: E=ZeV
*/
SECTION("High energy") {
// define projectile kinematics in lab frame
HEPMassType const projectileMass = 1_GeV;
HEPMomentumType P0 = 1_ZeV;
Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_PeV, -P0}};
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
const FourVector PprojLab(eProjectileLab, pProjectileLab);
// define boost to com frame
COMBoost boost(PprojLab, targetMass);
// boost projecticle
auto const PprojCoM = boost.toCoM(PprojLab);
// boost target
auto const PtargCoM = boost.toCoM(FourVector(targetMass, pTargetLab));
// sum of momenta in CoM, should be 0
auto const sumPCoM =
PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
CHECK(sumPCoM.norm() / P0 == Approx(0).margin(absMargin)); // MAKE RELATIVE CHECK
}
SECTION("rest frame") {
HEPMassType const projectileMass = 1_GeV;
HEPMomentumType const P0 = 1_TeV;
Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, P0, 0_GeV}};
HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
const FourVector PprojLab(eProjectileLab, pProjectileLab);
COMBoost boostRest(pProjectileLab, projectileMass);
auto const& csPrime = boostRest.GetRotatedCS();
FourVector const rest4Mom = boostRest.toCoM(PprojLab);
CHECK(rest4Mom.GetTimeLikeComponent() / 1_GeV == Approx(projectileMass / 1_GeV));
CHECK(rest4Mom.GetSpaceLikeComponents().norm() / 1_GeV ==
Approx(0).margin(absMargin));
FourVector const a{0_eV, Vector{csPrime, 0_eV, 5_GeV, 0_eV}};
FourVector const b{0_eV, Vector{rootCS, 3_GeV, 0_eV, 0_eV}};
auto const aLab = boostRest.fromCoM(a);
auto const bLab = boostRest.fromCoM(b);
CHECK(aLab.GetNorm() / a.GetNorm() == Approx(1));
CHECK(aLab.GetSpaceLikeComponents().GetComponents(csPrime)[1].magnitude() ==
Approx((5_GeV).magnitude()));
CHECK(bLab.GetSpaceLikeComponents().GetComponents(rootCS)[0].magnitude() ==
Approx((3_GeV).magnitude()));
}
}
...@@ -82,6 +82,7 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa ...@@ -82,6 +82,7 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa
plane_.GetNormal().dot(velocity.cross(magneticfield)) * 2 * k)) - plane_.GetNormal().dot(velocity.cross(magneticfield)) * 2 * k)) -
velocity.dot(plane_.GetNormal()) / velocity.GetNorm() ) / velocity.dot(plane_.GetNormal()) / velocity.GetNorm() ) /
(plane_.GetNormal().dot(velocity.cross(magneticfield)) * k); (plane_.GetNormal().dot(velocity.cross(magneticfield)) * k);
if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) { if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
return std::numeric_limits<double>::infinity() * 1_m; return std::numeric_limits<double>::infinity() * 1_m;
} else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) { } else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
......
/*
* (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/process/observation_plane/ObservationPlane.h>
#include <fstream>
using namespace corsika::process::observation_plane;
using namespace corsika::units::si;
ObservationPlane::ObservationPlane(
geometry::Plane const& obsPlane,
geometry::Vector<units::si::dimensionless_d> const& x_axis,
std::string const& filename, bool deleteOnHit)
: plane_(obsPlane)
, outputStream_(filename)
, deleteOnHit_(deleteOnHit)
, xAxis_(x_axis.normalized())
, yAxis_(obsPlane.GetNormal().cross(xAxis_)) {
outputStream_ << "#PDG code, energy / eV, x distance / m, y distance / m" << std::endl;
}
corsika::process::EProcessReturn ObservationPlane::DoContinuous(
setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) {
TimeType const timeOfIntersection =
(plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
trajectory.GetV0().dot(plane_.GetNormal());
if (timeOfIntersection < TimeType::zero()) { return process::EProcessReturn::eOk; }
if (plane_.IsAbove(trajectory.GetR0()) == plane_.IsAbove(trajectory.GetPosition(1))) {
return process::EProcessReturn::eOk;
}
auto const displacement = trajectory.GetPosition(1) - plane_.GetCenter();
outputStream_ << static_cast<int>(particles::GetPDG(particle.GetPID())) << ' '
<< particle.GetEnergy() * (1 / 1_eV) << ' '
<< displacement.dot(xAxis_) / 1_m << ' ' << displacement.dot(yAxis_) / 1_m
<< ' ' << std::endl;
if (deleteOnHit_) {
return process::EProcessReturn::eParticleAbsorbed;
} else {
return process::EProcessReturn::eOk;
}
}
LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
setup::Trajectory const& trajectory) {
TimeType const timeOfIntersection =
(plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
trajectory.GetV0().dot(plane_.GetNormal());
if (timeOfIntersection < TimeType::zero()) {
return std::numeric_limits<double>::infinity() * 1_m;
}
auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
}
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/FourVector.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h>
#include <corsika/process/sibyll/sibyll2.3d.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/utl/COMBoost.h>
#include <tuple>
using std::cout;
using std::endl;
using std::tuple;
using namespace corsika;
using namespace corsika::setup;
using SetupParticle = setup::Stack::StackIterator;
using SetupProjectile = setup::StackView::StackIterator;
using Track = Trajectory;
namespace corsika::process::sibyll {
Interaction::Interaction() {}
Interaction::~Interaction() {
cout << "Sibyll::Interaction n=" << count_ << " Nnuc=" << nucCount_ << endl;
}
void Interaction::Init() {
using random::RNGManager;
// initialize Sibyll
if (!initialized_) {
sibyll_ini_();
initialized_ = true;
}
}
void Interaction::SetAllStable() {
for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = -1 * abs(s_csydec_.idb[i]);
}
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
Interaction::GetCrossSection(const particles::Code BeamId,
const particles::Code TargetId,
const units::si::HEPEnergyType CoMenergy) const {
using namespace units::si;
double sigProd, sigEla, dummy, dum1, dum3, dum4;
double dumdif[3];
const int iBeam = process::sibyll::GetSibyllXSCode(BeamId);
if (!IsValidCoMEnergy(CoMenergy)) {
throw std::runtime_error(
"Interaction: GetCrossSection: CoM energy outside range for Sibyll!");
}
const double dEcm = CoMenergy / 1_GeV;
if (particles::IsNucleus(TargetId)) {
const int iTarget = particles::GetNucleusA(TargetId);
if (iTarget > maxTargetMassNumber_ || iTarget == 0)
throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 are allowed.");
sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
} else if (TargetId == particles::Proton::GetCode()) {
sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
} else {
// no interaction in sibyll possible, return infinite cross section? or throw?
sigProd = std::numeric_limits<double>::infinity();
sigEla = std::numeric_limits<double>::infinity();
}
return std::make_tuple(sigProd * 1_mb, sigEla * 1_mb);
}
template <>
units::si::GrammageType Interaction::GetInteractionLength(
SetupParticle const& vP) const {
using namespace units;
using namespace units::si;
using namespace geometry;
// coordinate system, get global frame of reference
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
const particles::Code corsikaBeamId = vP.GetPID();
// beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
const bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
// FOR NOW: assume target is at rest
MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
// total momentum and energy
HEPEnergyType Elab = vP.GetEnergy() + constants::nucleonMass;
MomentumVector pTotLab(rootCS, {0_GeV, 0_GeV, 0_GeV});
pTotLab += vP.GetMomentum();
pTotLab += pTarget;
auto const pTotLabNorm = pTotLab.norm();
// calculate cm. energy
const HEPEnergyType ECoM = sqrt(
(Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
cout << "Interaction: LambdaInt: \n"
<< " input energy: " << vP.GetEnergy() / 1_GeV << endl
<< " beam can interact:" << kInteraction << endl
<< " beam pid:" << vP.GetPID() << endl;
// TODO: move limits into variables
// FR: removed && Elab >= 8.5_GeV
if (kInteraction && IsValidCoMEnergy(ECoM)) {
// get target from environment
/*
the target should be defined by the Environment,
ideally as full particle object so that the four momenta
and the boosts can be defined..
*/
auto const* currentNode = vP.GetNode();
const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
si::CrossSectionType weightedProdCrossSection = mediumComposition.WeightedSum(
[=](particles::Code targetID) -> si::CrossSectionType {
return std::get<0>(this->GetCrossSection(corsikaBeamId, targetID, ECoM));
});
cout << "Interaction: "
<< "IntLength: weighted CrossSection (mb): " << weightedProdCrossSection / 1_mb
<< endl;
// calculate interaction length in medium
GrammageType const int_length = mediumComposition.GetAverageMassNumber() *
units::constants::u / weightedProdCrossSection;
cout << "Interaction: "
<< "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
<< endl;
return int_length;
}
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
/**
In this function SIBYLL is called to produce one event. The
event is copied (and boosted) into the shower lab frame.
*/
template <>
process::EProcessReturn Interaction::DoInteraction(SetupProjectile& vP) {
using namespace units;
using namespace utl;
using namespace units::si;
using namespace geometry;
const auto corsikaBeamId = vP.GetPID();
cout << "ProcessSibyll: "
<< "DoInteraction: " << corsikaBeamId << " interaction? "
<< process::sibyll::CanInteract(corsikaBeamId) << endl;
if (particles::IsNucleus(corsikaBeamId)) {
// nuclei handled by different process, this should not happen
throw std::runtime_error("Nuclear projectile are not handled by SIBYLL!");
}
if (process::sibyll::CanInteract(corsikaBeamId)) {
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// position and time of interaction, not used in Sibyll
Point pOrig = vP.GetPosition();
TimeType tOrig = vP.GetTime();
// define target
// for Sibyll is always a single nucleon
// FOR NOW: target is always at rest
const auto eTargetLab = 0_GeV + constants::nucleonMass;
const auto pTargetLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
const FourVector PtargLab(eTargetLab, pTargetLab);
// define projectile
HEPEnergyType const eProjectileLab = vP.GetEnergy();
auto const pProjectileLab = vP.GetMomentum();
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;
const FourVector PprojLab(eProjectileLab, pProjectileLab);
// define target kinematics in lab frame
// define boost to and from CoM frame
// CoM frame definition in Sibyll projectile: +z
COMBoost const boost(PprojLab, constants::nucleonMass);
// just for show:
// boost projecticle
auto const PprojCoM = boost.toCoM(PprojLab);
// boost target
auto const PtargCoM = boost.toCoM(PtargLab);
cout << "Interaction: ebeam CoM: " << PprojCoM.GetTimeLikeComponent() / 1_GeV
<< endl
<< "Interaction: pbeam CoM: "
<< PprojCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
cout << "Interaction: etarget CoM: " << PtargCoM.GetTimeLikeComponent() / 1_GeV
<< endl
<< "Interaction: ptarget CoM: "
<< PtargCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl;
cout << "Interaction: time: " << tOrig << endl;
HEPEnergyType Etot = eProjectileLab + eTargetLab;
MomentumVector Ptot = vP.GetMomentum();
// invariant mass, i.e. cm. energy
HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
// sample target mass number
auto const* currentNode = vP.GetNode();
auto const& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// get cross sections for target materials
/*
Here we read the cross section from the interaction model again,
should be passed from GetInteractionLength if possible
*/
//#warning reading interaction cross section again, should not be necessary
auto const& compVec = mediumComposition.GetComponents();
std::vector<si::CrossSectionType> cross_section_of_components(compVec.size());
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
const auto [sigProd, sigEla] = GetCrossSection(corsikaBeamId, targetId, Ecm);
[[maybe_unused]] const auto& dummy_sigEla = sigEla;
cross_section_of_components[i] = sigProd;
}
const auto targetCode =
mediumComposition.SampleTarget(cross_section_of_components, RNG_);
cout << "Interaction: target selected: " << targetCode << endl;
/*
FOR NOW: allow nuclei with A<18 or protons only.
when medium composition becomes more complex, approximations will have to be
allowed air in atmosphere also contains some Argon.
*/
int targetSibCode = -1;
if (IsNucleus(targetCode)) targetSibCode = GetNucleusA(targetCode);
if (targetCode == particles::Proton::GetCode()) targetSibCode = 1;
cout << "Interaction: sibyll code: " << targetSibCode << endl;
if (targetSibCode > maxTargetMassNumber_ || targetSibCode < 1)
throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 or protons are "
"allowed.");
// beam id for sibyll
const int kBeam = process::sibyll::ConvertToSibyllRaw(corsikaBeamId);
cout << "Interaction: "
<< " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
<< " Ecm(GeV): " << Ecm / 1_GeV << endl;
if (Ecm > GetMaxEnergyCoM())
throw std::runtime_error("Interaction::DoInteraction: CoM energy too high!");
// FR: removed eProjectileLab < 8.5_GeV ||
if (Ecm < GetMinEnergyCoM()) {
cout << "Interaction: "
<< " DoInteraction: should have dropped particle.. "
<< "THIS IS AN ERROR" << endl;
throw std::runtime_error("energy too low for SIBYLL");
} else {
count_++;
// Sibyll does not know about units..
const double sqs = Ecm / 1_GeV;
// running sibyll, filling stack
sibyll_(kBeam, targetSibCode, sqs);
// print final state
int print_unit = 6;
sib_list_(print_unit);
nucCount_ += get_nwounded() - 1;
// add particles from sibyll to stack
// link to sibyll stack
SibStack ss;
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
for (auto& psib : ss) {
// abort on particles that have decayed in Sibyll. Should not happen!
if (psib.HasDecayed())
throw std::runtime_error("found particle that decayed in SIBYLL!");
// transform energy to lab. frame
auto const pCoM = psib.GetMomentum();
HEPEnergyType const eCoM = psib.GetEnergy();
auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
// 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,
tOrig});
Plab_final += pnew.GetMomentum();
Elab_final += pnew.GetEnergy();
Ecm_final += psib.GetEnergy();