Skip to content

Commit 240aaff

Browse files
jesgumjesgum
committed
add particle production time to otf decayer
1 parent a8d638d commit 240aaff

5 files changed

Lines changed: 134 additions & 111 deletions

File tree

‎ALICE3/Core/Decayer.h‎

Lines changed: 18 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -61,28 +61,27 @@ class Decayer
6161
const double ctau = o2::constants::physics::LightSpeedCm2S * particleInfo->Lifetime(); // cm
6262
const double betaGamma = particle.p() / mass;
6363
const double rxyz = -betaGamma * ctau * std::log(1 - u);
64-
double vx, vy, vz;
6564
double px, py, e;
6665

6766
if (!charge) {
68-
vx = particle.vx() + rxyz * (particle.px() / particle.p());
69-
vy = particle.vy() + rxyz * (particle.py() / particle.p());
70-
vz = particle.vz() + rxyz * (particle.pz() / particle.p());
67+
mVx = particle.vx() + rxyz * (particle.px() / particle.p());
68+
mVy = particle.vy() + rxyz * (particle.py() / particle.p());
69+
mVz = particle.vz() + rxyz * (particle.pz() / particle.p());
7170
px = particle.px();
7271
py = particle.py();
7372
} else {
7473
o2::track::TrackParCov track;
7574
o2::math_utils::CircleXYf_t circle;
7675
o2::upgrade::convertOTFParticleToO2Track(particle, track, pdgDB);
7776

78-
float sna, csa;
77+
float sna{}, csa{};
7978
track.getCircleParams(mBz, circle, sna, csa);
8079
const double rxy = rxyz / std::sqrt(1. + track.getTgl() * track.getTgl());
8180
const double theta = rxy / circle.rC;
8281

83-
vx = ((particle.vx() - circle.xC) * std::cos(theta) - (particle.vy() - circle.yC) * std::sin(theta)) + circle.xC;
84-
vy = ((particle.vy() - circle.yC) * std::cos(theta) + (particle.vx() - circle.xC) * std::sin(theta)) + circle.yC;
85-
vz = particle.vz() + rxyz * (particle.pz() / track.getP());
82+
mVx = ((particle.vx() - circle.xC) * std::cos(theta) - (particle.vy() - circle.yC) * std::sin(theta)) + circle.xC;
83+
mVy = ((particle.vy() - circle.yC) * std::cos(theta) + (particle.vx() - circle.xC) * std::sin(theta)) + circle.yC;
84+
mVz = particle.vz() + rxyz * (particle.pz() / track.getP());
8685

8786
px = particle.px() * std::cos(theta) - particle.py() * std::sin(theta);
8887
py = particle.py() * std::cos(theta) + particle.px() * std::sin(theta);
@@ -125,25 +124,32 @@ class Decayer
125124
o2::upgrade::OTFParticle particle;
126125
TLorentzVector dau = *decay.GetDecay(i);
127126
particle.setPDG(pdgCodesDaughters[i]);
128-
particle.setVxVyVz(vx, vy, vz);
127+
particle.setVxVyVz(mVx, mVy, mVz);
129128
particle.setPxPyPzE(dau.Px(), dau.Py(), dau.Pz(), dau.E());
130129
decayProducts.push_back(particle);
131130
}
132131

133132
return decayProducts;
134133
}
135134

135+
// Setters
136+
void setBField(const double b) { mBz = b; }
136137
void setSeed(const int seed)
137138
{
138139
mRand3.SetSeed(seed); // For decay length sampling
139140
gRandom->SetSeed(seed); // For TGenPhaseSpace
140141
}
141142

142-
void setBField(const double b) { mBz = b; }
143+
// Getters
144+
float getSecondaryVertexX() const { return static_cast<float>(mVx); }
145+
float getSecondaryVertexY() const { return static_cast<float>(mVy); }
146+
float getSecondaryVertexZ() const { return static_cast<float>(mVz); }
147+
float getDecayRadius() const { return static_cast<float>(std::hypot(mVx, mVy)); }
143148

144149
private:
145-
TRandom3 mRand3;
146-
double mBz;
150+
double mBz{20.}; // kG
151+
double mVx{-1.}, mVy{-1.}, mVz{-1.};
152+
TRandom3 mRand3{};
147153
};
148154

149155
} // namespace upgrade

‎ALICE3/Core/TrackUtilities.cxx‎

Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -45,3 +45,71 @@ void o2::upgrade::convertTLorentzVectorToO2Track(const int charge,
4545
// Initialize TrackParCov in-place
4646
new (&o2track)(o2::track::TrackParCov)(x, particle.Phi(), params, covm);
4747
}
48+
49+
float o2::upgrade::computeParticleVelocity(float momentum, float mass)
50+
{
51+
const float a = momentum / mass;
52+
// uses light speed in cm/ps so output is in those units
53+
return o2::constants::physics::LightSpeedCm2PS * a / std::sqrt((1.f + a * a));
54+
}
55+
56+
float o2::upgrade::computeTrackLength(o2::track::TrackParCov track, float radius, float magneticField)
57+
{
58+
// don't make use of the track parametrization
59+
float length = -100;
60+
61+
o2::math_utils::CircleXYf_t trcCircle;
62+
float sna, csa;
63+
track.getCircleParams(magneticField, trcCircle, sna, csa);
64+
65+
// distance between circle centers (one circle is at origin -> easy)
66+
const float centerDistance = std::hypot(trcCircle.xC, trcCircle.yC);
67+
68+
// condition of circles touching - if not satisfied returned length will be -100
69+
if (centerDistance < trcCircle.rC + radius && centerDistance > std::fabs(trcCircle.rC - radius)) {
70+
length = 0.0f;
71+
72+
// base radical direction
73+
const float ux = trcCircle.xC / centerDistance;
74+
const float uy = trcCircle.yC / centerDistance;
75+
// calculate perpendicular vector (normalized) for +/- displacement
76+
const float vx = -uy;
77+
const float vy = +ux;
78+
// calculate coordinate for radical line
79+
const float radical = (centerDistance * centerDistance - trcCircle.rC * trcCircle.rC + radius * radius) / (2.0f * centerDistance);
80+
// calculate absolute displacement from center-to-center axis
81+
const float displace = (0.5f / centerDistance) * std::sqrt(
82+
(-centerDistance + trcCircle.rC - radius) *
83+
(-centerDistance - trcCircle.rC + radius) *
84+
(-centerDistance + trcCircle.rC + radius) *
85+
(centerDistance + trcCircle.rC + radius));
86+
87+
// possible intercept points of track and TOF layer in 2D plane
88+
const float point1[2] = {radical * ux + displace * vx, radical * uy + displace * vy};
89+
const float point2[2] = {radical * ux - displace * vx, radical * uy - displace * vy};
90+
91+
// decide on correct intercept point
92+
std::array<float, 3> mom;
93+
track.getPxPyPzGlo(mom);
94+
const float scalarProduct1 = point1[0] * mom[0] + point1[1] * mom[1];
95+
const float scalarProduct2 = point2[0] * mom[0] + point2[1] * mom[1];
96+
97+
// get start point
98+
std::array<float, 3> startPoint;
99+
track.getXYZGlo(startPoint);
100+
101+
float cosAngle = -1000, modulus = -1000;
102+
103+
if (scalarProduct1 > scalarProduct2) {
104+
modulus = std::hypot(point1[0] - trcCircle.xC, point1[1] - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC);
105+
cosAngle = (point1[0] - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (point1[1] - trcCircle.yC) * (startPoint[1] - trcCircle.yC);
106+
} else {
107+
modulus = std::hypot(point2[0] - trcCircle.xC, point2[1] - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC);
108+
cosAngle = (point2[0] - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (point2[1] - trcCircle.yC) * (startPoint[1] - trcCircle.yC);
109+
}
110+
cosAngle /= modulus;
111+
length = trcCircle.rC * std::acos(cosAngle);
112+
length *= std::sqrt(1.0f + track.getTgl() * track.getTgl());
113+
}
114+
return length;
115+
}

‎ALICE3/Core/TrackUtilities.h‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -104,6 +104,17 @@ o2::track::TrackParCov convertMCParticleToO2Track(McParticleType& particle,
104104
return o2track;
105105
}
106106

107+
/// returns velocity in centimeters per picoseconds
108+
/// \param momentum the momentum of the track
109+
/// \param mass the mass of the particle
110+
float computeParticleVelocity(float momentum, float mass);
111+
112+
/// function to calculate track length of this track up to a certain radius
113+
/// \param track the input track (TrackParCov)
114+
/// \param radius the radius of the layer you're calculating the length to
115+
/// \param magneticField the magnetic field to use when propagating
116+
float computeTrackLength(o2::track::TrackParCov track, float radius, float magneticField);
117+
107118
} // namespace o2::upgrade
108119

109120
#endif // ALICE3_CORE_TRACKUTILITIES_H_

‎ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx‎

Lines changed: 21 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,6 @@
1616
///
1717

1818
#include "ALICE3/Core/Decayer.h"
19-
#include "ALICE3/Core/OTFParticle.h"
2019
#include "ALICE3/Core/TrackUtilities.h"
2120
#include "ALICE3/DataModel/OTFMCParticle.h"
2221

@@ -77,16 +76,17 @@ struct OnTheFlyDecayer {
7776
o2::upgrade::Decayer decayer;
7877
Service<o2::framework::O2DatabasePDG> pdgDB;
7978
std::map<int, std::vector<o2::upgrade::OTFParticle>> mDecayDaughters;
79+
HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
8080

8181
Configurable<int> seed{"seed", 0, "Set seed for particle decayer"};
8282
Configurable<float> magneticField{"magneticField", 20., "Magnetic field (kG)"};
8383
Configurable<LabeledArray<int>> enabledDecays{"enabledDecays",
8484
{DefaultParameters[0], NumDecays, NumParameters, ParticleNames, ParameterNames},
8585
"Enable option for particle to be decayed: 0 - no, 1 - yes"};
8686

87-
HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
88-
87+
static constexpr float PicoToNano = 1.e-3f;
8988
int mCollisionId{-1};
89+
9090
std::vector<int> mEnabledDecays;
9191
void init(o2::framework::InitContext&)
9292
{
@@ -133,12 +133,29 @@ struct OnTheFlyDecayer {
133133

134134
particle.setIsAlive(false);
135135
std::vector<o2::upgrade::OTFParticle> decayStack = decayer.decayParticle(pdgDB, particle);
136+
const float decayRadius = decayer.getDecayRadius();
137+
const float trackVelocity = o2::upgrade::computeParticleVelocity(particle.p(), pdgDB->GetParticle(particle.pdgCode())->Mass());
138+
const int charge = pdgDB->GetParticle(particle.pdgCode())->Charge() / 3;
139+
float trackLength{-1.f};
140+
if (!charge) {
141+
const float dx = particle.vx() - decayer.getSecondaryVertexX();
142+
const float dy = particle.vy() - decayer.getSecondaryVertexY();
143+
const float dz = particle.vz() - decayer.getSecondaryVertexZ();
144+
trackLength = std::hypot(dx, dy, dz);
145+
} else {
146+
o2::track::TrackParCov o2track;
147+
o2::upgrade::convertOTFParticleToO2Track(particle, o2track, pdgDB);
148+
trackLength = o2::upgrade::computeTrackLength(o2track, decayRadius, magneticField);
149+
}
150+
151+
const float trackTimeNS = trackLength / trackVelocity * PicoToNano;
136152
particle.setIndicesDaughter(allParticles.size(), allParticles.size() + (decayStack.size() - 1));
137153
for (o2::upgrade::OTFParticle daughter : decayStack) {
138154
daughter.setIndicesMother(i, i);
139155
daughter.setCollisionId(mCollisionId);
140156
daughter.setIsAlive(true);
141157
daughter.setIsPrimary(false);
158+
daughter.setProductionTime(trackTimeNS);
142159
allParticles.push_back(daughter);
143160
ndau++;
144161
}
@@ -175,7 +192,7 @@ struct OnTheFlyDecayer {
175192
histos.fill(HIST("hNaNBookkeeping"), 0);
176193
}
177194

178-
// todo: status codes and vt
195+
// todo: status codes
179196
tableMcParticlesWithDau(otfParticle.collisionId(), otfParticle.pdgCode(), otfParticle.statusCode(),
180197
otfParticle.flags(), otfParticle.getMotherSpan(), otfParticle.getDaughters().data(), otfParticle.weight(),
181198
otfParticle.px(), otfParticle.py(), otfParticle.pz(), otfParticle.e(),

‎ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx‎

Lines changed: 16 additions & 95 deletions
Original file line numberDiff line numberDiff line change
@@ -448,81 +448,6 @@ struct OnTheFlyTofPid {
448448
return outerTOFLayer.isTrackInActiveArea(track);
449449
}
450450

451-
/// function to calculate track length of this track up to a certain radius
452-
/// \param track the input track
453-
/// \param radius the radius of the layer you're calculating the length to
454-
/// \param magneticField the magnetic field to use when propagating
455-
static float computeTrackLength(o2::track::TrackParCov track, float radius, float magneticField)
456-
{
457-
// don't make use of the track parametrization
458-
float length = -100;
459-
460-
o2::math_utils::CircleXYf_t trcCircle;
461-
float sna, csa;
462-
track.getCircleParams(magneticField, trcCircle, sna, csa);
463-
464-
// distance between circle centers (one circle is at origin -> easy)
465-
const float centerDistance = std::hypot(trcCircle.xC, trcCircle.yC);
466-
467-
// condition of circles touching - if not satisfied returned length will be -100
468-
if (centerDistance < trcCircle.rC + radius && centerDistance > std::fabs(trcCircle.rC - radius)) {
469-
length = 0.0f;
470-
471-
// base radical direction
472-
const float ux = trcCircle.xC / centerDistance;
473-
const float uy = trcCircle.yC / centerDistance;
474-
// calculate perpendicular vector (normalized) for +/- displacement
475-
const float vx = -uy;
476-
const float vy = +ux;
477-
// calculate coordinate for radical line
478-
const float radical = (centerDistance * centerDistance - trcCircle.rC * trcCircle.rC + radius * radius) / (2.0f * centerDistance);
479-
// calculate absolute displacement from center-to-center axis
480-
const float displace = (0.5f / centerDistance) * std::sqrt(
481-
(-centerDistance + trcCircle.rC - radius) *
482-
(-centerDistance - trcCircle.rC + radius) *
483-
(-centerDistance + trcCircle.rC + radius) *
484-
(centerDistance + trcCircle.rC + radius));
485-
486-
// possible intercept points of track and TOF layer in 2D plane
487-
const float point1[2] = {radical * ux + displace * vx, radical * uy + displace * vy};
488-
const float point2[2] = {radical * ux - displace * vx, radical * uy - displace * vy};
489-
490-
// decide on correct intercept point
491-
std::array<float, 3> mom;
492-
track.getPxPyPzGlo(mom);
493-
const float scalarProduct1 = point1[0] * mom[0] + point1[1] * mom[1];
494-
const float scalarProduct2 = point2[0] * mom[0] + point2[1] * mom[1];
495-
496-
// get start point
497-
std::array<float, 3> startPoint;
498-
track.getXYZGlo(startPoint);
499-
500-
float cosAngle = -1000, modulus = -1000;
501-
502-
if (scalarProduct1 > scalarProduct2) {
503-
modulus = std::hypot(point1[0] - trcCircle.xC, point1[1] - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC);
504-
cosAngle = (point1[0] - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (point1[1] - trcCircle.yC) * (startPoint[1] - trcCircle.yC);
505-
} else {
506-
modulus = std::hypot(point2[0] - trcCircle.xC, point2[1] - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC);
507-
cosAngle = (point2[0] - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (point2[1] - trcCircle.yC) * (startPoint[1] - trcCircle.yC);
508-
}
509-
cosAngle /= modulus;
510-
length = trcCircle.rC * std::acos(cosAngle);
511-
length *= std::sqrt(1.0f + track.getTgl() * track.getTgl());
512-
}
513-
return length;
514-
}
515-
516-
/// returns velocity in centimeters per picoseconds
517-
/// \param momentum the momentum of the tarck
518-
/// \param mass the mass of the particle
519-
float computeParticleVelocity(float momentum, float mass)
520-
{
521-
const float a = momentum / mass;
522-
// uses light speed in cm/ps so output is in those units
523-
return o2::constants::physics::LightSpeedCm2PS * a / std::sqrt((1.f + a * a));
524-
}
525-
526451
struct TracksWithTime {
527452
TracksWithTime(int pdgCode,
528453
std::pair<float, float> innerTOFTime,
@@ -696,8 +621,8 @@ struct OnTheFlyTofPid {
696621
}
697622
float trackLengthInnerTOF = -1, trackLengthOuterTOF = -1;
698623
if (xPv > kTrkXThreshold) {
699-
trackLengthInnerTOF = computeTrackLength(o2track, simConfig.innerTOFRadius, mMagneticField);
700-
trackLengthOuterTOF = computeTrackLength(o2track, simConfig.outerTOFRadius, mMagneticField);
624+
trackLengthInnerTOF = o2::upgrade::computeTrackLength(o2track, simConfig.innerTOFRadius, mMagneticField);
625+
trackLengthOuterTOF = o2::upgrade::computeTrackLength(o2track, simConfig.outerTOFRadius, mMagneticField);
701626
}
702627

703628
// Check if the track hit a sensitive area of the TOF
@@ -707,9 +632,7 @@ struct OnTheFlyTofPid {
707632
} else {
708633
float x = 0.f;
709634
o2track.getXatLabR(simConfig.innerTOFRadius, x, mMagneticField);
710-
if (plotsConfig.doQAplots) {
711-
histos.fill(HIST("iTOF/h2HitMap"), o2track.getZAt(x, mMagneticField), simConfig.innerTOFRadius * o2track.getPhiAt(x, mMagneticField));
712-
}
635+
histos.fill(HIST("iTOF/h2HitMap"), o2track.getZAt(x, mMagneticField), simConfig.innerTOFRadius * o2track.getPhiAt(x, mMagneticField));
713636
}
714637

715638
const bool activeOuterTOF = isInOuterTOFActiveArea(o2track);
@@ -718,9 +641,7 @@ struct OnTheFlyTofPid {
718641
} else {
719642
float x = 0.f;
720643
o2track.getXatLabR(simConfig.outerTOFRadius, x, mMagneticField);
721-
if (plotsConfig.doQAplots) {
722-
histos.fill(HIST("oTOF/h2HitMap"), o2track.getZAt(x, mMagneticField), simConfig.outerTOFRadius * o2track.getPhiAt(x, mMagneticField));
723-
}
644+
histos.fill(HIST("oTOF/h2HitMap"), o2track.getZAt(x, mMagneticField), simConfig.outerTOFRadius * o2track.getPhiAt(x, mMagneticField));
724645
}
725646

726647
// get mass to calculate velocity
@@ -730,7 +651,7 @@ struct OnTheFlyTofPid {
730651
upgradeTofMC(-999.f, -999.f, -999.f, -999.f);
731652
continue;
732653
}
733-
const float v = computeParticleVelocity(o2track.getP(), pdgInfo->Mass());
654+
const float v = o2::upgrade::computeParticleVelocity(o2track.getP(), pdgInfo->Mass());
734655
const float expectedTimeInnerTOF = trackLengthInnerTOF > 0 ? trackLengthInnerTOF / v + eventCollisionTimePS : -999.f; // arrival time to the Inner TOF in ps
735656
const float expectedTimeOuterTOF = trackLengthOuterTOF > 0 ? trackLengthOuterTOF / v + eventCollisionTimePS : -999.f; // arrival time to the Outer TOF in ps
736657
upgradeTofMC(expectedTimeInnerTOF, trackLengthInnerTOF, expectedTimeOuterTOF, trackLengthOuterTOF);
@@ -748,8 +669,8 @@ struct OnTheFlyTofPid {
748669
xPv = recoTrack.getX();
749670
}
750671
if (xPv > kTrkXThreshold) {
751-
trackLengthRecoInnerTOF = computeTrackLength(recoTrack, simConfig.innerTOFRadius, mMagneticField);
752-
trackLengthRecoOuterTOF = computeTrackLength(recoTrack, simConfig.outerTOFRadius, mMagneticField);
672+
trackLengthRecoInnerTOF = o2::upgrade::computeTrackLength(recoTrack, simConfig.innerTOFRadius, mMagneticField);
673+
trackLengthRecoOuterTOF = o2::upgrade::computeTrackLength(recoTrack, simConfig.outerTOFRadius, mMagneticField);
753674
}
754675

755676
// cache the track info needed for the event time calculation
@@ -870,7 +791,7 @@ struct OnTheFlyTofPid {
870791
nSigmaOuterTOF[ii] = -100;
871792

872793
momentumHypotheses[ii] = rigidity * kParticleCharges[ii]; // Total momentum for this hypothesis
873-
const float v = computeParticleVelocity(momentumHypotheses[ii], kParticleMasses[ii]);
794+
const float v = o2::upgrade::computeParticleVelocity(momentumHypotheses[ii], kParticleMasses[ii]);
874795

875796
expectedTimeInnerTOF[ii] = trackLengthInnerTOF / v;
876797
expectedTimeOuterTOF[ii] = trackLengthOuterTOF / v;
@@ -942,15 +863,15 @@ struct OnTheFlyTofPid {
942863
}
943864
}
944865
}
866+
}
945867

946-
const float deltaTrackLengthInnerTOF = std::abs(trackLengthInnerTOF - trackLengthRecoInnerTOF);
947-
if (trackLengthInnerTOF > 0 && trackLengthRecoInnerTOF > 0) {
948-
histos.fill(HIST("iTOF/h2dDeltaTrackLengthInnerVsPt"), noSmearingPt, deltaTrackLengthInnerTOF);
949-
}
950-
const float deltaTrackLengthOuterTOF = std::abs(trackLengthOuterTOF - trackLengthRecoOuterTOF);
951-
if (trackLengthOuterTOF > 0 && trackLengthRecoOuterTOF > 0) {
952-
histos.fill(HIST("oTOF/h2dDeltaTrackLengthOuterVsPt"), noSmearingPt, deltaTrackLengthOuterTOF);
953-
}
868+
const float deltaTrackLengthInnerTOF = std::abs(trackLengthInnerTOF - trackLengthRecoInnerTOF);
869+
if (trackLengthInnerTOF > 0 && trackLengthRecoInnerTOF > 0) {
870+
histos.fill(HIST("iTOF/h2dDeltaTrackLengthInnerVsPt"), noSmearingPt, deltaTrackLengthInnerTOF);
871+
}
872+
const float deltaTrackLengthOuterTOF = std::abs(trackLengthOuterTOF - trackLengthRecoOuterTOF);
873+
if (trackLengthOuterTOF > 0 && trackLengthRecoOuterTOF > 0) {
874+
histos.fill(HIST("oTOF/h2dDeltaTrackLengthOuterVsPt"), noSmearingPt, deltaTrackLengthOuterTOF);
954875
}
955876

956877
// Sigmas have been fully calculated. Please populate the NSigma helper table (once per track)

0 commit comments

Comments
 (0)