Skip to content

Commit 2180ed2

Browse files
committed
add production time to daughters in otf-decayer
1 parent 350dfb6 commit 2180ed2

6 files changed

Lines changed: 292 additions & 235 deletions

File tree

ALICE3/Core/Decayer.h

Lines changed: 19 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@
1919
#ifndef ALICE3_CORE_DECAYER_H_
2020
#define ALICE3_CORE_DECAYER_H_
2121

22+
#include "ALICE3/Core/OTFParticle.h"
2223
#include "ALICE3/Core/TrackUtilities.h"
2324

2425
#include <CommonConstants/PhysicsConstants.h>
@@ -60,28 +61,27 @@ class Decayer
6061
const double ctau = o2::constants::physics::LightSpeedCm2S * particleInfo->Lifetime(); // cm
6162
const double betaGamma = particle.p() / mass;
6263
const double rxyz = -betaGamma * ctau * std::log(1 - u);
63-
double vx, vy, vz;
6464
double px, py, e;
6565

6666
if (!charge) {
67-
vx = particle.vx() + rxyz * (particle.px() / particle.p());
68-
vy = particle.vy() + rxyz * (particle.py() / particle.p());
69-
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());
7070
px = particle.px();
7171
py = particle.py();
7272
} else {
7373
o2::track::TrackParCov track;
7474
o2::math_utils::CircleXYf_t circle;
7575
o2::upgrade::convertOTFParticleToO2Track(particle, track, pdgDB);
7676

77-
float sna, csa;
77+
float sna{}, csa{};
7878
track.getCircleParams(mBz, circle, sna, csa);
7979
const double rxy = rxyz / std::sqrt(1. + track.getTgl() * track.getTgl());
8080
const double theta = rxy / circle.rC;
8181

82-
vx = ((particle.vx() - circle.xC) * std::cos(theta) - (particle.vy() - circle.yC) * std::sin(theta)) + circle.xC;
83-
vy = ((particle.vy() - circle.yC) * std::cos(theta) + (particle.vx() - circle.xC) * std::sin(theta)) + circle.yC;
84-
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());
8585

8686
px = particle.px() * std::cos(theta) - particle.py() * std::sin(theta);
8787
py = particle.py() * std::cos(theta) + particle.px() * std::sin(theta);
@@ -124,25 +124,32 @@ class Decayer
124124
o2::upgrade::OTFParticle particle;
125125
TLorentzVector dau = *decay.GetDecay(i);
126126
particle.setPDG(pdgCodesDaughters[i]);
127-
particle.setVxVyVz(vx, vy, vz);
127+
particle.setVxVyVz(mVx, mVy, mVz);
128128
particle.setPxPyPzE(dau.Px(), dau.Py(), dau.Pz(), dau.E());
129129
decayProducts.push_back(particle);
130130
}
131131

132132
return decayProducts;
133133
}
134134

135+
// Setters
136+
void setBField(const double b) { mBz = b; }
135137
void setSeed(const int seed)
136138
{
137139
mRand3.SetSeed(seed); // For decay length sampling
138140
gRandom->SetSeed(seed); // For TGenPhaseSpace
139141
}
140142

141-
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)); }
142148

143149
private:
144-
TRandom3 mRand3;
145-
double mBz;
150+
double mBz{20.}; // kG
151+
double mVx{-1.}, mVy{-1.}, mVz{-1.};
152+
TRandom3 mRand3{};
146153
};
147154

148155
} // namespace upgrade

ALICE3/Core/OTFParticle.h

Lines changed: 165 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,165 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
12+
///
13+
/// \file OTFParticle.h
14+
/// \brief Basic class to hold information regarding a mc particle to be used in fast simulation
15+
/// \author Jesper Karlsson Gumprecht <jesper.gumprecht@cern.ch>
16+
///
17+
18+
#ifndef ALICE3_CORE_OTFPARTICLE_H_
19+
#define ALICE3_CORE_OTFPARTICLE_H_
20+
21+
#include <CommonConstants/MathConstants.h>
22+
23+
#include <array>
24+
#include <cmath>
25+
#include <cstdint>
26+
#include <span>
27+
28+
namespace o2::upgrade
29+
{
30+
31+
class OTFParticle
32+
{
33+
public:
34+
OTFParticle() = default;
35+
36+
template <typename TParticle>
37+
explicit OTFParticle(const TParticle& particle)
38+
{
39+
mPdgCode = particle.pdgCode();
40+
mGlobalIndex = particle.globalIndex();
41+
mCollisionId = particle.mcCollisionId();
42+
mPx = particle.px();
43+
mPy = particle.py();
44+
mPz = particle.pz();
45+
mE = particle.e();
46+
mVx = particle.vx();
47+
mVy = particle.vy();
48+
mVz = particle.vz();
49+
mVt = particle.vt();
50+
mIsFromMcParticles = true;
51+
if (particle.has_mothers()) {
52+
mIndicesMother = {particle.mothersIds().front(), particle.mothersIds().back()};
53+
}
54+
}
55+
56+
// Setters
57+
void setIsAlive(const bool isAlive) { mIsAlive = isAlive; }
58+
void setIsPrimary(const bool isPrimary) { mIsPrimary = isPrimary; }
59+
void setCollisionId(const int collisionId) { mCollisionId = collisionId; }
60+
void setPDG(const int pdg) { mPdgCode = pdg; }
61+
void setIndicesMother(const int start, const int stop) { mIndicesMother = {start, stop}; }
62+
void setIndicesDaughter(const int start, const int stop) { mIndicesDaughter = {start, stop}; }
63+
void setProductionTime(const float vt) { mVt = vt; }
64+
void setVxVyVz(const float vx, const float vy, const float vz)
65+
{
66+
mVx = vx;
67+
mVy = vy;
68+
mVz = vz;
69+
}
70+
void setPxPyPzE(const float px, const float py, const float pz, const float e)
71+
{
72+
mPx = px;
73+
mPy = py;
74+
mPz = pz;
75+
mE = e;
76+
}
77+
78+
// Getters
79+
int pdgCode() const { return mPdgCode; }
80+
int globalIndex() const { return mGlobalIndex; }
81+
int collisionId() const { return mCollisionId; }
82+
bool isAlive() const { return mIsAlive; }
83+
bool isPrimary() const { return mIsPrimary; }
84+
bool isFromMcParticles() const { return mIsFromMcParticles; }
85+
float weight() const
86+
{
87+
static constexpr float Weight = 1.f;
88+
return Weight;
89+
}
90+
uint8_t flags() const
91+
{
92+
static constexpr uint8_t Flags = 1;
93+
return Flags; // todo
94+
}
95+
int statusCode() const
96+
{
97+
static constexpr int StatusCode = 1;
98+
return StatusCode; // todo
99+
}
100+
float vx() const { return mVx; }
101+
float vy() const { return mVy; }
102+
float vz() const { return mVz; }
103+
float vt() const { return mVt; }
104+
float px() const { return mPx; }
105+
float py() const { return mPy; }
106+
float pz() const { return mPz; }
107+
float e() const { return mE; }
108+
float radius() const { return std::hypot(mVx, mVy); }
109+
float r() const { return radius(); }
110+
float pt() const { return std::hypot(mPx, mPy); }
111+
float p() const { return std::hypot(mPx, mPy, mPz); }
112+
float phi() const { return o2::constants::math::PI + std::atan2(-1.0f * py(), -1.0f * px()); }
113+
float eta() const
114+
{
115+
// As https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1943
116+
static constexpr float Tolerance = 1e-7f;
117+
if ((p() - mPz) < Tolerance) {
118+
return (mPz < 0.0f) ? -100.0f : 100.0f;
119+
} else {
120+
return 0.5f * std::log((p() + mPz) / (p() - mPz));
121+
}
122+
}
123+
float y() const
124+
{
125+
// As https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1922
126+
static constexpr float Tolerance = 1e-7f;
127+
if ((e() - mPz) < Tolerance) {
128+
return (mPz < 0.0f) ? -100.0f : 100.0f;
129+
} else {
130+
return 0.5f * std::log((mE + mPz) / (mE - mPz));
131+
}
132+
}
133+
int getMotherIndexStart() const { return mIndicesMother[0]; }
134+
int getMotherIndexStop() const { return mIndicesMother[1]; }
135+
int getDaughterIndexStart() const { return mIndicesDaughter[0]; }
136+
int getDaughterIndexStop() const { return mIndicesDaughter[1]; }
137+
std::array<int, 2> getMothers() const { return mIndicesMother; }
138+
std::array<int, 2> getDaughters() const { return mIndicesDaughter; }
139+
std::span<const int> getMotherSpan() const { return hasMothers() ? std::span<const int>(mIndicesMother.data(), 2) : std::span<const int>(); }
140+
141+
// Checks
142+
bool hasDaughters() const { return (mIndicesDaughter[0] > 0); }
143+
bool hasMothers() const { return (mIndicesMother[0] > 0); }
144+
bool hasNaN() const
145+
{
146+
return std::isnan(mPx) || std::isnan(mPy) || std::isnan(mPz) || std::isnan(mE) ||
147+
std::isnan(mVx) || std::isnan(mVy) || std::isnan(mVz);
148+
}
149+
bool hasIndex() const
150+
{
151+
return (mGlobalIndex != -1);
152+
}
153+
154+
private:
155+
int mPdgCode{}, mGlobalIndex{-1};
156+
int mCollisionId{};
157+
float mVx{}, mVy{}, mVz{}, mVt{};
158+
float mPx{}, mPy{}, mPz{}, mE{};
159+
bool mIsAlive{}, mIsFromMcParticles{false};
160+
bool mIsPrimary{};
161+
std::array<int, 2> mIndicesMother{-1, -1}, mIndicesDaughter{-1, -1};
162+
};
163+
164+
} // namespace o2::upgrade
165+
#endif // ALICE3_CORE_OTFPARTICLE_H_

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+
}

0 commit comments

Comments
 (0)