diff --git a/ALICE3/Core/Decayer.h b/ALICE3/Core/Decayer.h index 7189223573c..51296ffe67d 100644 --- a/ALICE3/Core/Decayer.h +++ b/ALICE3/Core/Decayer.h @@ -61,13 +61,12 @@ class Decayer const double ctau = o2::constants::physics::LightSpeedCm2S * particleInfo->Lifetime(); // cm const double betaGamma = particle.p() / mass; const double rxyz = -betaGamma * ctau * std::log(1 - u); - double vx, vy, vz; double px, py, e; if (!charge) { - vx = particle.vx() + rxyz * (particle.px() / particle.p()); - vy = particle.vy() + rxyz * (particle.py() / particle.p()); - vz = particle.vz() + rxyz * (particle.pz() / particle.p()); + mVx = particle.vx() + rxyz * (particle.px() / particle.p()); + mVy = particle.vy() + rxyz * (particle.py() / particle.p()); + mVz = particle.vz() + rxyz * (particle.pz() / particle.p()); px = particle.px(); py = particle.py(); } else { @@ -75,14 +74,14 @@ class Decayer o2::math_utils::CircleXYf_t circle; o2::upgrade::convertOTFParticleToO2Track(particle, track, pdgDB); - float sna, csa; + float sna{}, csa{}; track.getCircleParams(mBz, circle, sna, csa); const double rxy = rxyz / std::sqrt(1. + track.getTgl() * track.getTgl()); const double theta = rxy / circle.rC; - vx = ((particle.vx() - circle.xC) * std::cos(theta) - (particle.vy() - circle.yC) * std::sin(theta)) + circle.xC; - vy = ((particle.vy() - circle.yC) * std::cos(theta) + (particle.vx() - circle.xC) * std::sin(theta)) + circle.yC; - vz = particle.vz() + rxyz * (particle.pz() / track.getP()); + mVx = ((particle.vx() - circle.xC) * std::cos(theta) - (particle.vy() - circle.yC) * std::sin(theta)) + circle.xC; + mVy = ((particle.vy() - circle.yC) * std::cos(theta) + (particle.vx() - circle.xC) * std::sin(theta)) + circle.yC; + mVz = particle.vz() + rxyz * (particle.pz() / track.getP()); px = particle.px() * std::cos(theta) - particle.py() * std::sin(theta); py = particle.py() * std::cos(theta) + particle.px() * std::sin(theta); @@ -125,7 +124,7 @@ class Decayer o2::upgrade::OTFParticle particle; TLorentzVector dau = *decay.GetDecay(i); particle.setPDG(pdgCodesDaughters[i]); - particle.setVxVyVz(vx, vy, vz); + particle.setVxVyVz(mVx, mVy, mVz); particle.setPxPyPzE(dau.Px(), dau.Py(), dau.Pz(), dau.E()); decayProducts.push_back(particle); } @@ -133,17 +132,24 @@ class Decayer return decayProducts; } + // Setters + void setBField(const double b) { mBz = b; } void setSeed(const int seed) { mRand3.SetSeed(seed); // For decay length sampling gRandom->SetSeed(seed); // For TGenPhaseSpace } - void setBField(const double b) { mBz = b; } + // Getters + float getSecondaryVertexX() const { return static_cast(mVx); } + float getSecondaryVertexY() const { return static_cast(mVy); } + float getSecondaryVertexZ() const { return static_cast(mVz); } + float getDecayRadius() const { return static_cast(std::hypot(mVx, mVy)); } private: - TRandom3 mRand3; - double mBz; + double mBz{20.}; // kG + double mVx{-1.}, mVy{-1.}, mVz{-1.}; + TRandom3 mRand3{}; }; } // namespace upgrade diff --git a/ALICE3/Core/TrackUtilities.cxx b/ALICE3/Core/TrackUtilities.cxx index 739e1df28ae..4239864e933 100644 --- a/ALICE3/Core/TrackUtilities.cxx +++ b/ALICE3/Core/TrackUtilities.cxx @@ -45,3 +45,71 @@ void o2::upgrade::convertTLorentzVectorToO2Track(const int charge, // Initialize TrackParCov in-place new (&o2track)(o2::track::TrackParCov)(x, particle.Phi(), params, covm); } + +float o2::upgrade::computeParticleVelocity(float momentum, float mass) +{ + const float a = momentum / mass; + // uses light speed in cm/ps so output is in those units + return o2::constants::physics::LightSpeedCm2PS * a / std::sqrt((1.f + a * a)); +} + +float o2::upgrade::computeTrackLength(o2::track::TrackParCov track, float radius, float magneticField) +{ + // don't make use of the track parametrization + float length = -100; + + o2::math_utils::CircleXYf_t trcCircle; + float sna, csa; + track.getCircleParams(magneticField, trcCircle, sna, csa); + + // distance between circle centers (one circle is at origin -> easy) + const float centerDistance = std::hypot(trcCircle.xC, trcCircle.yC); + + // condition of circles touching - if not satisfied returned length will be -100 + if (centerDistance < trcCircle.rC + radius && centerDistance > std::fabs(trcCircle.rC - radius)) { + length = 0.0f; + + // base radical direction + const float ux = trcCircle.xC / centerDistance; + const float uy = trcCircle.yC / centerDistance; + // calculate perpendicular vector (normalized) for +/- displacement + const float vx = -uy; + const float vy = +ux; + // calculate coordinate for radical line + const float radical = (centerDistance * centerDistance - trcCircle.rC * trcCircle.rC + radius * radius) / (2.0f * centerDistance); + // calculate absolute displacement from center-to-center axis + const float displace = (0.5f / centerDistance) * std::sqrt( + (-centerDistance + trcCircle.rC - radius) * + (-centerDistance - trcCircle.rC + radius) * + (-centerDistance + trcCircle.rC + radius) * + (centerDistance + trcCircle.rC + radius)); + + // possible intercept points of track and TOF layer in 2D plane + const float point1[2] = {radical * ux + displace * vx, radical * uy + displace * vy}; + const float point2[2] = {radical * ux - displace * vx, radical * uy - displace * vy}; + + // decide on correct intercept point + std::array mom; + track.getPxPyPzGlo(mom); + const float scalarProduct1 = point1[0] * mom[0] + point1[1] * mom[1]; + const float scalarProduct2 = point2[0] * mom[0] + point2[1] * mom[1]; + + // get start point + std::array startPoint; + track.getXYZGlo(startPoint); + + float cosAngle = -1000, modulus = -1000; + + if (scalarProduct1 > scalarProduct2) { + modulus = std::hypot(point1[0] - trcCircle.xC, point1[1] - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC); + cosAngle = (point1[0] - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (point1[1] - trcCircle.yC) * (startPoint[1] - trcCircle.yC); + } else { + modulus = std::hypot(point2[0] - trcCircle.xC, point2[1] - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC); + cosAngle = (point2[0] - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (point2[1] - trcCircle.yC) * (startPoint[1] - trcCircle.yC); + } + cosAngle /= modulus; + length = trcCircle.rC * std::acos(cosAngle); + length *= std::sqrt(1.0f + track.getTgl() * track.getTgl()); + } + return length; +} diff --git a/ALICE3/Core/TrackUtilities.h b/ALICE3/Core/TrackUtilities.h index d61c72956cb..b80ea1b2341 100644 --- a/ALICE3/Core/TrackUtilities.h +++ b/ALICE3/Core/TrackUtilities.h @@ -104,6 +104,17 @@ o2::track::TrackParCov convertMCParticleToO2Track(McParticleType& particle, return o2track; } +/// returns velocity in centimeters per picoseconds +/// \param momentum the momentum of the track +/// \param mass the mass of the particle +float computeParticleVelocity(float momentum, float mass); + +/// function to calculate track length of this track up to a certain radius +/// \param track the input track (TrackParCov) +/// \param radius the radius of the layer you're calculating the length to +/// \param magneticField the magnetic field to use when propagating +float computeTrackLength(o2::track::TrackParCov track, float radius, float magneticField); + } // namespace o2::upgrade #endif // ALICE3_CORE_TRACKUTILITIES_H_ diff --git a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx index 38c51110feb..4bf984db338 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx @@ -77,6 +77,7 @@ struct OnTheFlyDecayer { o2::upgrade::Decayer decayer; Service pdgDB; std::map> mDecayDaughters; + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; Configurable seed{"seed", 0, "Set seed for particle decayer"}; Configurable magneticField{"magneticField", 20., "Magnetic field (kG)"}; @@ -84,9 +85,9 @@ struct OnTheFlyDecayer { {DefaultParameters[0], NumDecays, NumParameters, ParticleNames, ParameterNames}, "Enable option for particle to be decayed: 0 - no, 1 - yes"}; - HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - + static constexpr float PicoToNano = 1.e-3f; int mCollisionId{-1}; + std::vector mEnabledDecays; void init(o2::framework::InitContext&) { @@ -133,12 +134,29 @@ struct OnTheFlyDecayer { particle.setIsAlive(false); std::vector decayStack = decayer.decayParticle(pdgDB, particle); + const float decayRadius = decayer.getDecayRadius(); + const float trackVelocity = o2::upgrade::computeParticleVelocity(particle.p(), pdgDB->GetParticle(particle.pdgCode())->Mass()); + const int charge = pdgDB->GetParticle(particle.pdgCode())->Charge() / 3; + float trackLength{-1.f}; + if (!charge) { + const float dx = particle.vx() - decayer.getSecondaryVertexX(); + const float dy = particle.vy() - decayer.getSecondaryVertexY(); + const float dz = particle.vz() - decayer.getSecondaryVertexZ(); + trackLength = std::hypot(dx, dy, dz); + } else { + o2::track::TrackParCov o2track; + o2::upgrade::convertOTFParticleToO2Track(particle, o2track, pdgDB); + trackLength = o2::upgrade::computeTrackLength(o2track, decayRadius, magneticField); + } + + const float trackTimeNS = trackLength / trackVelocity * PicoToNano; particle.setIndicesDaughter(allParticles.size(), allParticles.size() + (decayStack.size() - 1)); for (o2::upgrade::OTFParticle daughter : decayStack) { daughter.setIndicesMother(i, i); daughter.setCollisionId(mCollisionId); daughter.setIsAlive(true); daughter.setIsPrimary(false); + daughter.setProductionTime(trackTimeNS); allParticles.push_back(daughter); ndau++; } @@ -175,7 +193,7 @@ struct OnTheFlyDecayer { histos.fill(HIST("hNaNBookkeeping"), 0); } - // todo: status codes and vt + // todo: status codes tableMcParticlesWithDau(otfParticle.collisionId(), otfParticle.pdgCode(), otfParticle.statusCode(), otfParticle.flags(), otfParticle.getMotherSpan(), otfParticle.getDaughters().data(), otfParticle.weight(), otfParticle.px(), otfParticle.py(), otfParticle.pz(), otfParticle.e(), diff --git a/ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx index 4617ca73657..99d9de97f22 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx @@ -448,81 +448,6 @@ struct OnTheFlyTofPid { return outerTOFLayer.isTrackInActiveArea(track); } - /// function to calculate track length of this track up to a certain radius - /// \param track the input track - /// \param radius the radius of the layer you're calculating the length to - /// \param magneticField the magnetic field to use when propagating - static float computeTrackLength(o2::track::TrackParCov track, float radius, float magneticField) - { - // don't make use of the track parametrization - float length = -100; - - o2::math_utils::CircleXYf_t trcCircle; - float sna, csa; - track.getCircleParams(magneticField, trcCircle, sna, csa); - - // distance between circle centers (one circle is at origin -> easy) - const float centerDistance = std::hypot(trcCircle.xC, trcCircle.yC); - - // condition of circles touching - if not satisfied returned length will be -100 - if (centerDistance < trcCircle.rC + radius && centerDistance > std::fabs(trcCircle.rC - radius)) { - length = 0.0f; - - // base radical direction - const float ux = trcCircle.xC / centerDistance; - const float uy = trcCircle.yC / centerDistance; - // calculate perpendicular vector (normalized) for +/- displacement - const float vx = -uy; - const float vy = +ux; - // calculate coordinate for radical line - const float radical = (centerDistance * centerDistance - trcCircle.rC * trcCircle.rC + radius * radius) / (2.0f * centerDistance); - // calculate absolute displacement from center-to-center axis - const float displace = (0.5f / centerDistance) * std::sqrt( - (-centerDistance + trcCircle.rC - radius) * - (-centerDistance - trcCircle.rC + radius) * - (-centerDistance + trcCircle.rC + radius) * - (centerDistance + trcCircle.rC + radius)); - - // possible intercept points of track and TOF layer in 2D plane - const float point1[2] = {radical * ux + displace * vx, radical * uy + displace * vy}; - const float point2[2] = {radical * ux - displace * vx, radical * uy - displace * vy}; - - // decide on correct intercept point - std::array mom; - track.getPxPyPzGlo(mom); - const float scalarProduct1 = point1[0] * mom[0] + point1[1] * mom[1]; - const float scalarProduct2 = point2[0] * mom[0] + point2[1] * mom[1]; - - // get start point - std::array startPoint; - track.getXYZGlo(startPoint); - - float cosAngle = -1000, modulus = -1000; - - if (scalarProduct1 > scalarProduct2) { - modulus = std::hypot(point1[0] - trcCircle.xC, point1[1] - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC); - cosAngle = (point1[0] - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (point1[1] - trcCircle.yC) * (startPoint[1] - trcCircle.yC); - } else { - modulus = std::hypot(point2[0] - trcCircle.xC, point2[1] - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC); - cosAngle = (point2[0] - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (point2[1] - trcCircle.yC) * (startPoint[1] - trcCircle.yC); - } - cosAngle /= modulus; - length = trcCircle.rC * std::acos(cosAngle); - length *= std::sqrt(1.0f + track.getTgl() * track.getTgl()); - } - return length; - } - - /// returns velocity in centimeters per picoseconds - /// \param momentum the momentum of the tarck - /// \param mass the mass of the particle - float computeParticleVelocity(float momentum, float mass) - { - const float a = momentum / mass; - // uses light speed in cm/ps so output is in those units - return o2::constants::physics::LightSpeedCm2PS * a / std::sqrt((1.f + a * a)); - } - struct TracksWithTime { TracksWithTime(int pdgCode, std::pair innerTOFTime, @@ -696,8 +621,8 @@ struct OnTheFlyTofPid { } float trackLengthInnerTOF = -1, trackLengthOuterTOF = -1; if (xPv > kTrkXThreshold) { - trackLengthInnerTOF = computeTrackLength(o2track, simConfig.innerTOFRadius, mMagneticField); - trackLengthOuterTOF = computeTrackLength(o2track, simConfig.outerTOFRadius, mMagneticField); + trackLengthInnerTOF = o2::upgrade::computeTrackLength(o2track, simConfig.innerTOFRadius, mMagneticField); + trackLengthOuterTOF = o2::upgrade::computeTrackLength(o2track, simConfig.outerTOFRadius, mMagneticField); } // Check if the track hit a sensitive area of the TOF @@ -730,7 +655,7 @@ struct OnTheFlyTofPid { upgradeTofMC(-999.f, -999.f, -999.f, -999.f); continue; } - const float v = computeParticleVelocity(o2track.getP(), pdgInfo->Mass()); + const float v = o2::upgrade::computeParticleVelocity(o2track.getP(), pdgInfo->Mass()); const float expectedTimeInnerTOF = trackLengthInnerTOF > 0 ? trackLengthInnerTOF / v + eventCollisionTimePS : -999.f; // arrival time to the Inner TOF in ps const float expectedTimeOuterTOF = trackLengthOuterTOF > 0 ? trackLengthOuterTOF / v + eventCollisionTimePS : -999.f; // arrival time to the Outer TOF in ps upgradeTofMC(expectedTimeInnerTOF, trackLengthInnerTOF, expectedTimeOuterTOF, trackLengthOuterTOF); @@ -748,8 +673,8 @@ struct OnTheFlyTofPid { xPv = recoTrack.getX(); } if (xPv > kTrkXThreshold) { - trackLengthRecoInnerTOF = computeTrackLength(recoTrack, simConfig.innerTOFRadius, mMagneticField); - trackLengthRecoOuterTOF = computeTrackLength(recoTrack, simConfig.outerTOFRadius, mMagneticField); + trackLengthRecoInnerTOF = o2::upgrade::computeTrackLength(recoTrack, simConfig.innerTOFRadius, mMagneticField); + trackLengthRecoOuterTOF = o2::upgrade::computeTrackLength(recoTrack, simConfig.outerTOFRadius, mMagneticField); } // cache the track info needed for the event time calculation @@ -870,7 +795,7 @@ struct OnTheFlyTofPid { nSigmaOuterTOF[ii] = -100; momentumHypotheses[ii] = rigidity * kParticleCharges[ii]; // Total momentum for this hypothesis - const float v = computeParticleVelocity(momentumHypotheses[ii], kParticleMasses[ii]); + const float v = o2::upgrade::computeParticleVelocity(momentumHypotheses[ii], kParticleMasses[ii]); expectedTimeInnerTOF[ii] = trackLengthInnerTOF / v; expectedTimeOuterTOF[ii] = trackLengthOuterTOF / v;