Skip to content

Commit 483988d

Browse files
jesgumjesgum
committed
Add smearing to shortlived particles in LUT process function in otf tracker
1 parent 86ccf3c commit 483988d

3 files changed

Lines changed: 74 additions & 46 deletions

File tree

‎ALICE3/Core/Decayer.h‎

Lines changed: 54 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -31,13 +31,12 @@
3131
#include <TLorentzVector.h>
3232
#include <TRandom3.h>
3333

34+
#include <array>
3435
#include <cmath>
3536
#include <cstddef>
3637
#include <vector>
3738

38-
namespace o2
39-
{
40-
namespace upgrade
39+
namespace o2::upgrade
4140
{
4241

4342
class Decayer
@@ -47,44 +46,27 @@ class Decayer
4746
Decayer() = default;
4847

4948
template <typename TDatabase>
50-
std::vector<o2::upgrade::OTFParticle> decayParticle(const TDatabase& pdgDB, const OTFParticle& particle)
49+
std::vector<o2::upgrade::OTFParticle> decayParticle(const OTFParticle& particle, const TDatabase& pdgDB)
5150
{
52-
const auto& particleInfo = pdgDB->GetParticle(particle.pdgCode());
51+
auto particleInfo = pdgDB->GetParticle(particle.pdgCode());
5352
if (!particleInfo) {
5453
return {};
5554
}
5655

5756
const int charge = particleInfo->Charge() / 3;
5857
const double mass = particleInfo->Mass();
59-
60-
const double u = mRand3.Uniform(0.001, 0.999);
61-
const double ctau = o2::constants::physics::LightSpeedCm2S * particleInfo->Lifetime(); // cm
62-
const double betaGamma = particle.p() / mass;
63-
const double rxyz = -betaGamma * ctau * std::log(1 - u);
64-
double px, py, e;
58+
std::array<double, 3> decayVtx = generateDecayVertex<double>(particle, pdgDB);
59+
mVx = decayVtx[0];
60+
mVy = decayVtx[1];
61+
mVz = decayVtx[2];
62+
double px{}, py{}, e{};
6563

6664
if (!charge) {
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());
7065
px = particle.px();
7166
py = particle.py();
7267
} else {
73-
o2::track::TrackParCov track;
74-
o2::math_utils::CircleXYf_t circle;
75-
o2::upgrade::convertOTFParticleToO2Track(particle, track, pdgDB);
76-
77-
float sna{}, csa{};
78-
track.getCircleParams(mBz, circle, sna, csa);
79-
const double rxy = rxyz / std::sqrt(1. + track.getTgl() * track.getTgl());
80-
const double theta = rxy / circle.rC;
81-
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());
85-
86-
px = particle.px() * std::cos(theta) - particle.py() * std::sin(theta);
87-
py = particle.py() * std::cos(theta) + particle.px() * std::sin(theta);
68+
px = particle.px() * std::cos(mTheta) - particle.py() * std::sin(mTheta);
69+
py = particle.py() * std::cos(mTheta) + particle.px() * std::sin(mTheta);
8870
}
8971

9072
double brTotal = 0.;
@@ -133,6 +115,42 @@ class Decayer
133115
return decayProducts;
134116
}
135117

118+
template <typename T = float, typename TDatabase, typename TParticle>
119+
std::array<T, 3> generateDecayVertex(const TParticle& particle, const TDatabase& pdgDB)
120+
{
121+
std::array<T, 3> decayVertex{};
122+
auto particleInfo = pdgDB->GetParticle(particle.pdgCode());
123+
if (!particleInfo) {
124+
return {};
125+
}
126+
127+
const int charge = particleInfo->Charge() / 3;
128+
const double mass = particleInfo->Mass();
129+
const double u = mRand3.Uniform(0.001, 0.999);
130+
const double ctau = o2::constants::physics::LightSpeedCm2S * particleInfo->Lifetime(); // cm
131+
const double betaGamma = particle.p() / mass;
132+
const double rxyz = -betaGamma * ctau * std::log(1 - u);
133+
134+
if (!charge) {
135+
decayVertex[0] = particle.vx() + rxyz * (particle.px() / particle.p());
136+
decayVertex[1] = particle.vy() + rxyz * (particle.py() / particle.p());
137+
decayVertex[2] = particle.vz() + rxyz * (particle.pz() / particle.p());
138+
} else {
139+
o2::math_utils::CircleXYf_t circle;
140+
o2::track::TrackParCov track = o2::upgrade::convertMCParticleToO2Track(particle, pdgDB);
141+
142+
float sna{}, csa{};
143+
track.getCircleParams(mBz, circle, sna, csa);
144+
const double rxy = rxyz / std::sqrt(1. + track.getTgl() * track.getTgl());
145+
mTheta = rxy / circle.rC;
146+
147+
decayVertex[0] = ((particle.vx() - circle.xC) * std::cos(mTheta) - (particle.vy() - circle.yC) * std::sin(mTheta)) + circle.xC;
148+
decayVertex[1] = ((particle.vy() - circle.yC) * std::cos(mTheta) + (particle.vx() - circle.xC) * std::sin(mTheta)) + circle.yC;
149+
decayVertex[2] = particle.vz() + rxyz * (particle.pz() / track.getP());
150+
}
151+
return decayVertex;
152+
}
153+
136154
// Setters
137155
void setBField(const double b) { mBz = b; }
138156
void setSeed(const int seed)
@@ -142,18 +160,18 @@ class Decayer
142160
}
143161

144162
// 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)); }
163+
[[nodiscard]] float getSecondaryVertexX() const { return static_cast<float>(mVx); }
164+
[[nodiscard]] float getSecondaryVertexY() const { return static_cast<float>(mVy); }
165+
[[nodiscard]] float getSecondaryVertexZ() const { return static_cast<float>(mVz); }
166+
[[nodiscard]] float getDecayRadius() const { return static_cast<float>(std::hypot(mVx, mVy)); }
149167

150168
private:
151169
double mBz{20.}; // kG
152170
double mVx{-1.}, mVy{-1.}, mVz{-1.};
153-
TRandom3 mRand3{};
171+
double mTheta{};
172+
TRandom3 mRand3;
154173
};
155174

156-
} // namespace upgrade
157-
} // namespace o2
175+
} // namespace o2::upgrade
158176

159177
#endif // ALICE3_CORE_DECAYER_H_

‎ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -152,7 +152,7 @@ struct OnTheFlyDecayer {
152152
}
153153

154154
particle.setBitOff(o2::upgrade::DecayerBits::IsAlive);
155-
std::vector<o2::upgrade::OTFParticle> decayStack = decayer.decayParticle(pdgDB, particle);
155+
std::vector<o2::upgrade::OTFParticle> decayStack = decayer.decayParticle(particle, pdgDB);
156156
if (decayStack.empty()) {
157157
continue;
158158
}

‎ALICE3/TableProducer/OTF/onTheFlyTracker.cxx‎

Lines changed: 19 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,7 @@
2323
/// \author Roberto Preghenella preghenella@bo.infn.it
2424
///
2525

26+
#include "ALICE3/Core/Decayer.h"
2627
#include "ALICE3/Core/DetLayer.h"
2728
#include "ALICE3/Core/FastTracker.h"
2829
#include "ALICE3/Core/FlatTrackSmearer.h"
@@ -378,6 +379,9 @@ struct OnTheFlyTracker {
378379

379380
// Track smearer array, one per geometry
380381
std::vector<std::unique_ptr<o2::delphes::TrackSmearer>> mSmearer;
382+
// Configuration defined at init time
383+
o2::fastsim::GeometryContainer mGeoContainer;
384+
float mMagneticField = 0.0f;
381385

382386
// For processing and vertexing
383387
std::vector<TrackAlice3> recoPrimaries;
@@ -395,10 +399,8 @@ struct OnTheFlyTracker {
395399
// For TGenPhaseSpace seed
396400
TRandom3 rand;
397401
Service<o2::ccdb::BasicCCDBManager> ccdb{};
402+
o2::upgrade::Decayer decayer;
398403

399-
// Configuration defined at init time
400-
o2::fastsim::GeometryContainer mGeoContainer;
401-
float mMagneticField = 0.0f;
402404
// Time resolution constants
403405
static constexpr float timeResolutionNs = 100.f; // ns
404406
static constexpr float nsToMus = 1e-3f;
@@ -438,6 +440,7 @@ struct OnTheFlyTracker {
438440

439441
const int nGeometries = mGeoContainer.getNumberOfConfigurations();
440442
mMagneticField = mGeoContainer.getFloatValue(0, "global", "magneticfield");
443+
decayer.setBField(mMagneticField);
441444
for (int icfg = 0; icfg < nGeometries; ++icfg) {
442445
const std::string histPath = "Configuration_" + std::to_string(icfg) + "/";
443446
mSmearer.emplace_back(std::make_unique<o2::delphes::TrackSmearer>());
@@ -1910,7 +1913,6 @@ struct OnTheFlyTracker {
19101913
uint32_t multiplicityCounter = 0;
19111914
// Now that the multiplicity is known, we can process the particles to smear them
19121915
for (const auto& mcParticle : mcParticles) {
1913-
19141916
if (!mcParticle.isPhysicalPrimary()) {
19151917
continue;
19161918
}
@@ -1950,15 +1952,24 @@ struct OnTheFlyTracker {
19501952
bool reconstructed = true;
19511953
int nTrkHits = 0;
19521954
if (enablePrimarySmearing) {
1953-
if (fastPrimaryTrackerSettings.fastTrackPrimaries || fastPrimaryTrackerSettings.fastTrackShortLivedParticles) {
1954-
o2::track::TrackParCov perfectTrackParCov;
1955-
o2::upgrade::convertMCParticleToO2Track(mcParticle, perfectTrackParCov, pdgDB);
1955+
if (fastPrimaryTrackerSettings.fastTrackPrimaries && longLivedToBeHandled) {
1956+
o2::track::TrackParCov perfectTrackParCov = o2::upgrade::convertMCParticleToO2Track(mcParticle, pdgDB);
19561957
perfectTrackParCov.setPID(pdgCodeToPID(mcParticle.pdgCode()));
19571958
computeBremsstrahlungLoss(icfg, mcParticle, perfectTrackParCov);
19581959
nTrkHits = fastTracker[icfg]->FastTrack(perfectTrackParCov, trackParCov, dNdEta);
19591960
if (nTrkHits < fastPrimaryTrackerSettings.minSiliconHits) {
19601961
reconstructed = false;
19611962
}
1963+
} else if (fastPrimaryTrackerSettings.fastTrackShortLivedParticles && shortLivedToBeHandled) {
1964+
o2::track::TrackParCov perfectTrackParCov = o2::upgrade::convertMCParticleToO2Track(mcParticle, pdgDB);
1965+
perfectTrackParCov.setPID(pdgCodeToPID(mcParticle.pdgCode()));
1966+
computeBremsstrahlungLoss(icfg, mcParticle, perfectTrackParCov);
1967+
const std::array<float, 3> decayVtx = decayer.generateDecayVertex(mcParticle, pdgDB);
1968+
const float decayRadius2D = std::hypot(decayVtx[0], decayVtx[1]);
1969+
nTrkHits = fastTracker[icfg]->FastTrack(perfectTrackParCov, trackParCov, dNdEta, decayRadius2D);
1970+
if (nTrkHits < fastPrimaryTrackerSettings.minSiliconHits) {
1971+
reconstructed = false;
1972+
}
19621973
} else {
19631974
o2::upgrade::convertMCParticleToO2Track(mcParticle, trackParCov, pdgDB);
19641975
computeBremsstrahlungLoss(icfg, mcParticle, trackParCov);
@@ -2155,8 +2166,7 @@ struct OnTheFlyTracker {
21552166
computeBremsstrahlungLoss(icfg, mcParticle, trackParCov);
21562167
reconstructed = mSmearer[icfg]->smearTrack(trackParCov, mcParticle.pdgCode(), dNdEta);
21572168
} else if (shortLivedToBeHandled && fastPrimaryTrackerSettings.fastTrackShortLivedParticles) {
2158-
o2::track::TrackParCov perfectTrackParCov;
2159-
o2::upgrade::convertMCParticleToO2Track(mcParticle, perfectTrackParCov, pdgDB);
2169+
o2::track::TrackParCov perfectTrackParCov = o2::upgrade::convertMCParticleToO2Track(mcParticle, pdgDB);
21602170
perfectTrackParCov.setPID(pdgCodeToPID(mcParticle.pdgCode()));
21612171
computeBremsstrahlungLoss(icfg, mcParticle, perfectTrackParCov);
21622172
nTrkHits = fastTracker[icfg]->FastTrack(perfectTrackParCov, trackParCov, dNdEta, mcParticle.decayRadius());

0 commit comments

Comments
 (0)