Skip to content

Commit 53eaa79

Browse files
authored
Merge branch 'AliceO2Group:master' into master
2 parents e2acbfb + 2dd7aca commit 53eaa79

247 files changed

Lines changed: 30324 additions & 10094 deletions

File tree

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

.github/workflows/mega-linter.yml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -38,7 +38,7 @@ jobs:
3838
id: ml
3939
# You can override MegaLinter flavor used to have faster performances
4040
# More info at https://megalinter.io/flavors/
41-
uses: oxsecurity/megalinter@v9.4.0
41+
uses: oxsecurity/megalinter@v9.5.0
4242
env:
4343
# All available variables are described in documentation:
4444
# https://megalinter.io/configuration/

.mega-linter.yml

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,11 +21,13 @@ DISABLE_LINTERS:
2121
- REPOSITORY_DEVSKIM
2222
- REPOSITORY_GITLEAKS
2323
- REPOSITORY_KICS
24+
- REPOSITORY_OSV_SCANNER
2425
- REPOSITORY_SECRETLINT
2526
- REPOSITORY_TRIVY
2627
- YAML_PRETTIER
2728
- YAML_V8R
2829
DISABLE_ERRORS_LINTERS: # If errors are found by these linters, they will be considered as non blocking.
30+
- ACTION_ZIZMOR
2931
- PYTHON_BANDIT # The bandit check is overly broad and complains about subprocess usage.
3032
SHOW_ELAPSED_TIME: true
3133
FILEIO_REPORTER: false
@@ -42,3 +44,4 @@ CPP_CLANG_FORMAT_FILE_EXTENSIONS: [".C", ".c", ".c++", ".cc", ".cl", ".cpp", ".c
4244
CPP_CPPCHECK_FILE_EXTENSIONS: [".C", ".c", ".c++", ".cc", ".cl", ".cpp", ".cu", ".cuh", ".cxx", ".cxx.in", ".h", ".h++", ".hh", ".h.in", ".hpp", ".hxx", ".inc", ".inl", ".macro"]
4345
CPP_CPPCHECK_ARGUMENTS: --language=c++ --std=c++20 --check-level=exhaustive --suppressions-list=cppcheck_config
4446
REPOSITORY_GITLEAKS_PR_COMMITS_SCAN: true
47+
ACTION_ZIZMOR_UNSECURED_ENV_VARIABLES: [GITHUB_TOKEN]

ALICE3/Core/Decayer.h

Lines changed: 20 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,33 @@ 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());
129+
particle.setBitOn(o2::upgrade::DecayerBits::ProducedByDecayer);
129130
decayProducts.push_back(particle);
130131
}
131132

132133
return decayProducts;
133134
}
134135

136+
// Setters
137+
void setBField(const double b) { mBz = b; }
135138
void setSeed(const int seed)
136139
{
137140
mRand3.SetSeed(seed); // For decay length sampling
138141
gRandom->SetSeed(seed); // For TGenPhaseSpace
139142
}
140143

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

143150
private:
144-
TRandom3 mRand3;
145-
double mBz;
151+
double mBz{20.}; // kG
152+
double mVx{-1.}, mVy{-1.}, mVz{-1.};
153+
TRandom3 mRand3{};
146154
};
147155

148156
} // namespace upgrade

ALICE3/Core/OTFParticle.h

Lines changed: 187 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,187 @@
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 <bitset>
25+
#include <cmath>
26+
#include <cstdint>
27+
#include <span>
28+
29+
namespace o2::upgrade
30+
{
31+
32+
enum class DecayerBits { ProducedByDecayer = 0,
33+
IsPrimary,
34+
IsAlive };
35+
36+
class OTFParticle
37+
{
38+
public:
39+
OTFParticle() = default;
40+
41+
template <typename TParticle>
42+
explicit OTFParticle(const TParticle& particle)
43+
{
44+
mPdgCode = particle.pdgCode();
45+
mGlobalIndex = particle.globalIndex();
46+
mCollisionId = particle.mcCollisionId();
47+
mPx = particle.px();
48+
mPy = particle.py();
49+
mPz = particle.pz();
50+
mE = particle.e();
51+
mVx = particle.vx();
52+
mVy = particle.vy();
53+
mVz = particle.vz();
54+
mVt = particle.vt();
55+
mFlag = particle.flags();
56+
mStatusCode = particle.statusCode();
57+
mIsFromMcParticles = true;
58+
if (particle.has_mothers()) {
59+
mIndicesMother = {particle.mothersIds().front(), particle.mothersIds().back()};
60+
}
61+
if constexpr (requires { particle.decayerBits(); }) {
62+
mBits = particle.decayerBits();
63+
} else {
64+
// If we are here, we created particle in the standard workflow -- without secondaries
65+
// Then we should set all particles as physical primaries accordingly
66+
setBitOn(DecayerBits::IsPrimary);
67+
}
68+
}
69+
70+
// Setters
71+
void setIsPrimary(const bool isPrimary) { mIsPrimary = isPrimary; }
72+
void setCollisionId(const int collisionId) { mCollisionId = collisionId; }
73+
void setPDG(const int pdg) { mPdgCode = pdg; }
74+
void setIndicesMother(const int start, const int stop) { mIndicesMother = {start, stop}; }
75+
void setIndicesDaughter(const int start, const int stop) { mIndicesDaughter = {start, stop}; }
76+
void setProductionTime(const float vt) { mVt = vt; }
77+
void setFlags(uint8_t flag) { mFlag = flag; }
78+
void setVxVyVz(const float vx, const float vy, const float vz)
79+
{
80+
mVx = vx;
81+
mVy = vy;
82+
mVz = vz;
83+
}
84+
void setPxPyPzE(const float px, const float py, const float pz, const float e)
85+
{
86+
mPx = px;
87+
mPy = py;
88+
mPz = pz;
89+
mE = e;
90+
}
91+
92+
// Getters
93+
int pdgCode() const { return mPdgCode; }
94+
int globalIndex() const { return mGlobalIndex; }
95+
int collisionId() const { return mCollisionId; }
96+
bool isAlive() const { return mIsAlive; }
97+
bool isPrimary() const { return mIsPrimary; }
98+
bool isFromMcParticles() const { return mIsFromMcParticles; }
99+
float weight() const
100+
{
101+
static constexpr float Weight = 1.f;
102+
return Weight;
103+
}
104+
uint8_t flags() const { return mFlag; }
105+
int statusCode() const { return mStatusCode; }
106+
float vx() const { return mVx; }
107+
float vy() const { return mVy; }
108+
float vz() const { return mVz; }
109+
float vt() const { return mVt; }
110+
float px() const { return mPx; }
111+
float py() const { return mPy; }
112+
float pz() const { return mPz; }
113+
float e() const { return mE; }
114+
float radius() const { return std::hypot(mVx, mVy); }
115+
float r() const { return radius(); }
116+
float pt() const { return std::hypot(mPx, mPy); }
117+
float p() const { return std::hypot(mPx, mPy, mPz); }
118+
float phi() const { return o2::constants::math::PI + std::atan2(-1.0f * py(), -1.0f * px()); }
119+
float eta() const
120+
{
121+
// Conditionally defined to avoid FPEs
122+
// As https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1959
123+
static constexpr float Tolerance = 1e-7f;
124+
if ((p() - mPz) < Tolerance) {
125+
return (mPz < 0.0f) ? -100.0f : 100.0f;
126+
} else {
127+
return 0.5f * std::log((p() + mPz) / (p() - mPz));
128+
}
129+
}
130+
float y() const
131+
{
132+
// Conditionally defined to avoid FPEs
133+
// As https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1980
134+
static constexpr float Tolerance = 1e-7f;
135+
if ((e() - mPz) < Tolerance) {
136+
return (mPz < 0.0f) ? -100.0f : 100.0f;
137+
} else {
138+
return 0.5f * std::log((mE + mPz) / (mE - mPz));
139+
}
140+
}
141+
int getMotherIndexStart() const { return mIndicesMother[0]; }
142+
int getMotherIndexStop() const { return mIndicesMother[1]; }
143+
int getDaughterIndexStart() const { return mIndicesDaughter[0]; }
144+
int getDaughterIndexStop() const { return mIndicesDaughter[1]; }
145+
std::array<int, 2> getMothers() const { return mIndicesMother; }
146+
std::array<int, 2> getDaughters() const { return mIndicesDaughter; }
147+
std::span<const int> getMotherSpan() const { return hasMothers() ? std::span<const int>(mIndicesMother.data(), 2) : std::span<const int>(); }
148+
149+
// Checks
150+
bool hasDaughters() const { return (mIndicesDaughter[0] > 0); }
151+
bool hasMothers() const { return (mIndicesMother[0] > 0); }
152+
bool hasNaN() const
153+
{
154+
return std::isnan(mPx) || std::isnan(mPy) || std::isnan(mPz) || std::isnan(mE) ||
155+
std::isnan(mVx) || std::isnan(mVy) || std::isnan(mVz);
156+
}
157+
bool hasIndex() const
158+
{
159+
return (mGlobalIndex != -1);
160+
}
161+
162+
// Bits
163+
bool checkBit(DecayerBits bit) const { return mBits.test(static_cast<size_t>(bit)); }
164+
void setBit(DecayerBits bit, bool value = true) { mBits.set(static_cast<size_t>(bit), value); }
165+
void setBitOn(DecayerBits bit) { mBits.set(static_cast<size_t>(bit), true); }
166+
void setBitOff(DecayerBits bit) { mBits.set(static_cast<size_t>(bit), false); }
167+
168+
std::bitset<8> getBits() const { return mBits; }
169+
uint8_t getBitsValue() const { return static_cast<uint8_t>(mBits.to_ulong()); }
170+
void setBits(std::bitset<8> bits) { mBits = bits; }
171+
172+
private:
173+
int mPdgCode{}, mGlobalIndex{-1};
174+
int mCollisionId{};
175+
float mVx{}, mVy{}, mVz{}, mVt{};
176+
float mPx{}, mPy{}, mPz{}, mE{};
177+
bool mIsAlive{}, mIsFromMcParticles{false};
178+
bool mIsPrimary{};
179+
180+
int mStatusCode{};
181+
uint8_t mFlag{};
182+
std::bitset<8> mBits{};
183+
std::array<int, 2> mIndicesMother{-1, -1}, mIndicesDaughter{-1, -1};
184+
};
185+
186+
} // namespace o2::upgrade
187+
#endif // ALICE3_CORE_OTFPARTICLE_H_

ALICE3/Core/TrackUtilities.cxx

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,8 @@
1717

1818
#include "TrackUtilities.h"
1919

20+
#include <CommonConstants/PhysicsConstants.h>
21+
#include <MathUtils/Primitive2D.h>
2022
#include <MathUtils/Utils.h>
2123
#include <ReconstructionDataFormats/Track.h>
2224

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

0 commit comments

Comments
 (0)