Skip to content

Commit 8973025

Browse files
jesgumjesgum
committed
First step towards general handling of otf decays
1 parent 51db1b8 commit 8973025

File tree

3 files changed

+300
-1
lines changed

3 files changed

+300
-1
lines changed

‎ALICE3/Core/Decayer.h‎

Lines changed: 140 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,140 @@
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 Decayer.h
14+
/// \author Jesper Karlsson Gumprecht
15+
/// \since 15/12/2025
16+
///
17+
18+
#ifndef ALICE3_CORE_DECAYER_H_
19+
#define ALICE3_CORE_DECAYER_H_
20+
21+
#include "ALICE3/Core/TrackUtilities.h"
22+
#include "ReconstructionDataFormats/Track.h"
23+
24+
#include <TDatabasePDG.h>
25+
#include <TDecayChannel.h>
26+
#include <TGenPhaseSpace.h>
27+
#include <TLorentzVector.h>
28+
#include <TParticlePDG.h>
29+
#include <TRandom3.h>
30+
31+
#include <string>
32+
#include <unordered_map>
33+
#include <vector>
34+
#include <array>
35+
#include <cmath>
36+
37+
38+
namespace o2
39+
{
40+
namespace upgrade
41+
{
42+
43+
class Decayer {
44+
public:
45+
// Default constructor
46+
Decayer() = default;
47+
48+
template <typename TDatabase>
49+
std::vector<o2::upgrade::OTFParticle> decayParticle(const TDatabase& pdgDB, const o2::track::TrackParCov& track, const int pdgCode)
50+
{
51+
const auto& particleInfo = pdgDB->GetParticle(pdgCode);
52+
const int charge = particleInfo->Charge() / 3;
53+
const double mass = particleInfo->Mass();
54+
std::array<float, 3> mom;
55+
std::array<float, 3> pos;
56+
track.getPxPyPzGlo(mom);
57+
track.getXYZGlo(pos);
58+
59+
const double u = mRand3.Uniform(0, 1);
60+
const double ctau = o2::constants::physics::LightSpeedCm2S * particleInfo->Lifetime(); // cm
61+
const double betaGamma = track.getP() / mass;
62+
const double rxyz = -betaGamma * ctau * std::log(1 - u);
63+
64+
float sna, csa;
65+
o2::math_utils::CircleXYf_t circle;
66+
track.getCircleParams(mBz, circle, sna, csa);
67+
const double rxy = rxyz / std::sqrt(1. + track.getTgl() * track.getTgl());
68+
const double theta = rxy / circle.rC;
69+
70+
const double vx = ((pos[0] - circle.xC) * std::cos(theta) - (pos[1] - circle.yC) * std::sin(theta)) + circle.xC;
71+
const double vy = ((pos[1] - circle.yC) * std::cos(theta) + (pos[0] - circle.xC) * std::sin(theta)) + circle.yC;
72+
const double vz = mom[2] + rxyz * (mom[2] / track.getP());
73+
74+
const double px = mom[0] * std::cos(theta) - mom[1] * std::sin(theta);
75+
const double py = mom[1] * std::cos(theta) + mom[0] * std::sin(theta);
76+
const double e = std::sqrt(mass * mass + px * px + py * py + mom[2] * mom[2]);
77+
78+
double brTotal = 0.;
79+
for (int ch = 0; ch < particleInfo->NDecayChannels(); ++ch) {
80+
brTotal += particleInfo->DecayChannel(ch)->BranchingRatio();
81+
}
82+
83+
double brSum = 0.;
84+
std::vector<double> dauMasses;
85+
std::vector<int> pdgCodesDaughters;
86+
const double randomChannel = mRand3.Uniform(0., brTotal);
87+
for (int ch = 0; ch < particleInfo->NDecayChannels(); ++ch) {
88+
brSum += particleInfo->DecayChannel(ch)->BranchingRatio();
89+
if (randomChannel < brSum) {
90+
for (int dau = 0; dau < particleInfo->DecayChannel(ch)->NDaughters(); ++dau) {
91+
const int pdgDau = particleInfo->DecayChannel(ch)->DaughterPdgCode(dau);
92+
pdgCodesDaughters.push_back(pdgDau);
93+
const auto& dauInfo = pdgDB->GetParticle(pdgDau);
94+
dauMasses.push_back(dauInfo->Mass());
95+
}
96+
break;
97+
}
98+
}
99+
100+
TLorentzVector tlv(px, py, mom[2], e);
101+
TGenPhaseSpace decay;
102+
decay.SetDecay(tlv, dauMasses.size(), dauMasses.data());
103+
decay.Generate();
104+
105+
std::vector<o2::upgrade::OTFParticle> decayProducts;
106+
for (size_t i = 0; i < dauMasses.size(); ++i) {
107+
o2::upgrade::OTFParticle particle;
108+
TLorentzVector dau = *decay.GetDecay(i);
109+
particle.setPDG(pdgCodesDaughters[i]);
110+
particle.setVxVyVz(vx, vy, vz);
111+
particle.setPxPyPzE(dau.Px(), dau.Py(), dau.Pz(), dau.E());
112+
decayProducts.push_back(particle);
113+
}
114+
115+
return decayProducts;
116+
}
117+
118+
119+
// Setters
120+
void setSeed(const int seed)
121+
{
122+
mRand3.SetSeed(seed); // For decay length sampling
123+
gRandom->SetSeed(seed); // For TGenPhaseSpace
124+
}
125+
126+
void setBField(const double b) { mBz = b; }
127+
128+
// Getters
129+
130+
131+
132+
private:
133+
TRandom3 mRand3;
134+
double mBz;
135+
};
136+
137+
} // namespace fastsim
138+
} // namespace o2
139+
140+
#endif // ALICE3_CORE_DECAYER_H_

‎ALICE3/Core/TrackUtilities.h‎

Lines changed: 31 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -21,9 +21,10 @@
2121
#include "ReconstructionDataFormats/Track.h"
2222

2323
#include "TLorentzVector.h"
24-
2524
#include <vector>
2625

26+
27+
2728
namespace o2::upgrade
2829
{
2930

@@ -37,6 +38,17 @@ void convertTLorentzVectorToO2Track(const int charge,
3738
const std::vector<double> productionVertex,
3839
o2::track::TrackParCov& o2track);
3940

41+
struct OTFParticle {
42+
int pdgCode;
43+
float e;
44+
float vx, vy, vz;
45+
float px, py, pz;
46+
47+
void setPDG(int pdg) { pdgCode = pdg; }
48+
void setVxVyVz(float _vx, float _vy, float _vz) { vx = _vx; vy = _vy; vz = _vz; }
49+
void setPxPyPzE(float _px, float _py, float _pz, float _e) { px = _px; py = _py; pz = _pz; e = _e; }
50+
};
51+
4052
/// Function to convert a TLorentzVector into a perfect Track
4153
/// \param pdgCode particle pdg
4254
/// \param particle the particle to convert (TLorentzVector)
@@ -58,6 +70,24 @@ void convertTLorentzVectorToO2Track(int pdgCode,
5870
convertTLorentzVectorToO2Track(charge, particle, productionVertex, o2track);
5971
}
6072

73+
/// Function to convert a OTFParticle into a perfect Track
74+
/// \param particle the particle to convert (OTFParticle)
75+
/// \param o2track the address of the resulting TrackParCov
76+
/// \param pdg the pdg service
77+
template <typename PdgService>
78+
void convertOTFParticleToO2Track(const OTFParticle& particle, o2::track::TrackParCov& o2track, const PdgService& pdg)
79+
{
80+
const auto pdgInfo = pdg->GetParticle(particle.pdgCode);
81+
int charge = 0;
82+
if (pdgInfo != nullptr) {
83+
charge = pdgInfo->Charge() / 3;
84+
}
85+
86+
static TLorentzVector tlv;
87+
tlv.SetPxPyPzE(particle.px, particle.py, particle.pz, particle.e);
88+
convertTLorentzVectorToO2Track(particle.pdgCode, tlv, {particle.vx, particle.vy, particle.vz}, o2track, pdg);
89+
}
90+
6191
/// Function to convert a McParticle into a perfect Track
6292
/// \param particle the particle to convert (mcParticle)
6393
/// \param o2track the address of the resulting TrackParCov

‎ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx‎

Lines changed: 129 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,129 @@
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 onTheFlyDecayer.cxx
14+
/// \brief pre-processing for on-the-fly analysis
15+
/// \author Jesper Karlsson Gumprecht <jesper.gumprecht@cern.ch>
16+
///
17+
18+
#include "ALICE3/Core/Decayer.h"
19+
#include "ALICE3/Core/TrackUtilities.h"
20+
21+
#include <CCDB/BasicCCDBManager.h>
22+
#include <CommonConstants/MathConstants.h>
23+
#include <CommonConstants/PhysicsConstants.h>
24+
#include <DCAFitter/DCAFitterN.h>
25+
#include <DataFormatsParameters/GRPMagField.h>
26+
#include <DetectorsBase/Propagator.h>
27+
#include <DetectorsVertexing/PVertexer.h>
28+
#include <DetectorsVertexing/PVertexerHelpers.h>
29+
#include <Field/MagneticField.h>
30+
#include <Framework/AnalysisDataModel.h>
31+
#include <Framework/AnalysisTask.h>
32+
#include <Framework/HistogramRegistry.h>
33+
#include <Framework/O2DatabasePDGPlugin.h>
34+
#include <Framework/StaticFor.h>
35+
#include <Framework/runDataProcessing.h>
36+
#include <ReconstructionDataFormats/DCA.h>
37+
#include <SimulationDataFormat/InteractionSampler.h>
38+
39+
#include <TPDGCode.h>
40+
41+
#include <string>
42+
#include <vector>
43+
44+
using namespace o2;
45+
using namespace o2::framework;
46+
47+
static constexpr int kNumDecays = 7;
48+
static constexpr int kNumParameters = 1;
49+
static constexpr int defaultParameters[kNumDecays][kNumParameters]{{1}, {1}, {1}, {1}, {1}, {1}, {1}};
50+
static const std::vector<std::string> parameterNames{"enable"};
51+
static const std::vector<std::string> particleNames{"K0s",
52+
"Lambda",
53+
"Anti-Lambda",
54+
"Xi",
55+
"Anti-Xi",
56+
"Omega",
57+
"Anti-Omega"};
58+
static const std::vector<int> pdgCodes{kK0Short,
59+
kLambda0,
60+
kLambda0Bar,
61+
kXiMinus,
62+
kXiPlusBar,
63+
kOmegaMinus,
64+
kOmegaPlusBar};
65+
66+
struct OnTheFlyDecayer {
67+
o2::upgrade::Decayer decayer;
68+
Service<o2::ccdb::BasicCCDBManager> ccdb;
69+
Service<o2::framework::O2DatabasePDG> pdgDB;
70+
71+
Configurable<int> seed{"seed", 0, "Set seed for particle decayer"};
72+
Configurable<float> magneticField{"magneticField", 20., "Magnetic field (kG)"};
73+
Configurable<LabeledArray<int>> enabledDecays{"enabledDecays",
74+
{defaultParameters[0], kNumDecays, kNumParameters, particleNames, parameterNames},
75+
"Enable option for particle to be decayed: 0 - no, 1 - yes"};
76+
77+
std::vector<int> mEnabledDecays;
78+
void init(o2::framework::InitContext&)
79+
{
80+
ccdb->setURL("http://alice-ccdb.cern.ch");
81+
ccdb->setTimestamp(-1);
82+
decayer.setSeed(seed);
83+
decayer.setBField(magneticField);
84+
for (int i = 0; i < kNumDecays; ++i) {
85+
if (enabledDecays->get(particleNames[i].c_str(), "enable")) {
86+
mEnabledDecays.push_back(pdgCodes[i]);
87+
}
88+
}
89+
}
90+
91+
bool canDecay(const int pdgCode)
92+
{
93+
return std::find(mEnabledDecays.begin(), mEnabledDecays.end(), pdgCode) != mEnabledDecays.end();
94+
}
95+
96+
void process(aod::McCollision const&, aod::McParticles const& mcParticles)
97+
{
98+
for (const auto& particle : mcParticles) {
99+
if (!canDecay(particle.pdgCode())) {
100+
continue;
101+
}
102+
103+
o2::track::TrackParCov o2track;
104+
o2::upgrade::convertMCParticleToO2Track(particle, o2track, pdgDB);
105+
std::vector<o2::upgrade::OTFParticle> decayDaughters = decayer.decayParticle(pdgDB, o2track, particle.pdgCode());
106+
std::cout << particle.pdgCode() << " decayed into: ";
107+
for (const auto& dau : decayDaughters) {
108+
o2::track::TrackParCov dauTrack;
109+
o2::upgrade::convertOTFParticleToO2Track(dau, dauTrack, pdgDB);
110+
if (canDecay(dau.pdgCode)) {
111+
std::vector<o2::upgrade::OTFParticle> cascadingDaughers = decayer.decayParticle(pdgDB, dauTrack, dau.pdgCode);
112+
for (const auto& daudau : cascadingDaughers) {
113+
decayDaughters.push_back(daudau);
114+
}
115+
}
116+
}
117+
std::vector<bool> isFinalState;
118+
for (size_t i = 0; i < decayDaughters.size(); ++i) {
119+
isFinalState.push_back(!canDecay(decayDaughters[i].pdgCode));
120+
}
121+
}
122+
}
123+
};
124+
125+
126+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
127+
{
128+
return WorkflowSpec{adaptAnalysisTask<OnTheFlyDecayer>(cfgc)};
129+
}

0 commit comments

Comments
 (0)