Skip to content

Commit 4f42fb0

Browse files
jesgumjesgum
committed
add otf tracker process function to do new mcparticle table with daus
1 parent 398c722 commit 4f42fb0

File tree

5 files changed

+370
-46
lines changed

5 files changed

+370
-46
lines changed

‎ALICE3/Core/Decayer.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -69,7 +69,7 @@ class Decayer
6969
vy = pos[1] + rxyz * (mom[1] / track.getP());
7070
vz = pos[2] + rxyz * (mom[2] / track.getP());
7171
px = mom[0];
72-
py = mom[2];
72+
py = mom[1];
7373
} else {
7474
float sna, csa;
7575
o2::math_utils::CircleXYf_t circle;

‎ALICE3/Core/TrackUtilities.h‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -33,8 +33,10 @@ struct OTFParticle {
3333
float mE;
3434
float mVx, mVy, mVz;
3535
float mPx, mPy, mPz;
36+
bool mIsAlive;
3637

3738
// Setters
39+
void setIsAlive(bool isAlive) { mIsAlive = isAlive; }
3840
void setPDG(int pdg) { mPdgCode = pdg; }
3941
void setVxVyVz(float vx, float vy, float vz)
4042
{
@@ -52,6 +54,7 @@ struct OTFParticle {
5254

5355
// Getters
5456
int pdgCode() const { return mPdgCode; }
57+
int isAlive() const { return mIsAlive; }
5558
float vx() const { return mVx; }
5659
float vy() const { return mVy; }
5760
float vz() const { return mVz; }
@@ -97,7 +100,9 @@ void convertTLorentzVectorToO2Track(int pdgCode,
97100
/// \param o2track the address of the resulting TrackParCov
98101
/// \param pdg the pdg service
99102
template <typename PdgService>
100-
void convertOTFParticleToO2Track(const OTFParticle& particle, o2::track::TrackParCov& o2track, const PdgService& pdg)
103+
void convertOTFParticleToO2Track(const OTFParticle& particle,
104+
o2::track::TrackParCov& o2track,
105+
const PdgService& pdg)
101106
{
102107
static TLorentzVector tlv;
103108
tlv.SetPxPyPzE(particle.px(), particle.py(), particle.pz(), particle.e());

‎ALICE3/DataModel/OTFMCParticle.h‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,9 @@ DECLARE_SOA_COLUMN(Eta, eta, float);
3232
DECLARE_SOA_COLUMN(Pt, pt, float);
3333
DECLARE_SOA_COLUMN(P, p, float);
3434
DECLARE_SOA_COLUMN(Y, y, float);
35+
DECLARE_SOA_COLUMN(IsAlive, isAlive, bool);
36+
DECLARE_SOA_COLUMN(IsPrimary, isPrimary, bool);
37+
3538
} // namespace otfmcparticle
3639

3740
DECLARE_SOA_TABLE_FULL(McParticlesWithDau, "McParticlesWithDau", "AOD", "MCPARTICLEWITHDAU",
@@ -56,6 +59,8 @@ DECLARE_SOA_TABLE_FULL(McParticlesWithDau, "McParticlesWithDau", "AOD", "MCPARTI
5659
otfmcparticle::Pt,
5760
otfmcparticle::P,
5861
otfmcparticle::Y,
62+
otfmcparticle::IsAlive,
63+
otfmcparticle::IsPrimary,
5964
mcparticle::PVector<mcparticle::Px, mcparticle::Py, mcparticle::Pz>,
6065
mcparticle::ProducedByGenerator<mcparticle::Flags>,
6166
mcparticle::FromBackgroundEvent<mcparticle::Flags>,

‎ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx‎

Lines changed: 75 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -51,7 +51,7 @@ using namespace o2::framework;
5151
static constexpr int NumDecays = 7;
5252
static constexpr int NumParameters = 1;
5353
static constexpr int DefaultParameters[NumDecays][NumParameters]{{1}, {1}, {1}, {1}, {1}, {1}, {1}};
54-
static constexpr float VerySmall = 1e-7f;
54+
static constexpr float Tolerance = 1e-7f;
5555
static const std::vector<std::string> parameterNames{"enable"};
5656
static const std::vector<std::string> particleNames{"K0s",
5757
"Lambda",
@@ -102,6 +102,12 @@ struct OnTheFlyDecayer {
102102
histos.add("AntiXi/hGenAntiXi", "hGenAntiXi", kTH1D, {axisPt});
103103
histos.add("Omega/hGenOmega", "hGenOmega", kTH1D, {axisPt});
104104
histos.add("AntiOmega/hGenAntiOmega", "hGenAntiOmega", kTH1D, {axisPt});
105+
106+
histos.add("GeneratedElectron/hGenEl", "hGenEl", kTH1D, {axisPt});
107+
histos.add("GeneratedMuon/hGenMu", "hGenMu", kTH1D, {axisPt});
108+
histos.add("GeneratedPion/hGenPi", "hGenPi", kTH1D, {axisPt});
109+
histos.add("GeneratedKaon/hGenKa", "hGenKa", kTH1D, {axisPt});
110+
histos.add("GeneratedProton/hGenPr", "hGenPr", kTH1D, {axisPt});
105111
}
106112

107113
bool canDecay(const int pdgCode)
@@ -123,40 +129,61 @@ struct OnTheFlyDecayer {
123129
o2::upgrade::convertMCParticleToO2Track(particle, o2track, pdgDB);
124130
decayDaughters = decayer.decayParticle(pdgDB, o2track, particle.pdgCode());
125131
for (size_t idau{0}; idau < decayDaughters.size(); ++idau) {
126-
const auto& dau = decayDaughters[idau];
132+
o2::upgrade::OTFParticle& dau = decayDaughters[idau];
127133
o2::track::TrackParCov dauTrack;
128134
o2::upgrade::convertOTFParticleToO2Track(dau, dauTrack, pdgDB);
129135
if (canDecay(dau.pdgCode())) {
136+
dau.setIsAlive(false);
130137
std::vector<o2::upgrade::OTFParticle> cascadingDaughers = decayer.decayParticle(pdgDB, dauTrack, dau.pdgCode());
131138
for (const auto& daudau : cascadingDaughers) {
132139
decayDaughters.push_back(daudau);
133140
if (MaxNestedDecays < ++nDecays) {
134141
LOG(error) << "Seemingly stuck trying to perpetually decay products from pdg: " << particle.pdgCode();
135142
}
136143
}
144+
} else {
145+
dau.setIsAlive(true);
137146
}
138147
}
139148

140149
if (decayDaughters.empty()) {
141150
LOG(error) << "Attempted to decay " << particle.pdgCode() << " but resulting vector of daugthers were empty";
142151
continue;
143152
}
144-
}
145153

146-
if (particle.pdgCode() == kK0Short) {
147-
histos.fill(HIST("K0S/hGenK0S"), particle.pt());
148-
} else if (particle.pdgCode() == kLambda0) {
149-
histos.fill(HIST("Lambda/hGenLambda"), particle.pt());
150-
} else if (particle.pdgCode() == kLambda0Bar) {
151-
histos.fill(HIST("AntiLambda/hGenAntiLambda"), particle.pt());
152-
} else if (particle.pdgCode() == kXiMinus) {
153-
histos.fill(HIST("Xi/hGenXi"), particle.pt());
154-
} else if (particle.pdgCode() == kXiPlusBar) {
155-
histos.fill(HIST("AntiXi/hGenAntiXi"), particle.pt());
156-
} else if (particle.pdgCode() == kOmegaMinus) {
157-
histos.fill(HIST("Omega/hGenOmega"), particle.pt());
158-
} else if (particle.pdgCode() == kOmegaPlusBar) {
159-
histos.fill(HIST("AntiOmega/hGenAntiOmega"), particle.pt());
154+
155+
switch (particle.pdgCode()) {
156+
case kK0Short:
157+
histos.fill(HIST("K0S/hGenK0S"), particle.pt());
158+
break;
159+
160+
case kLambda0:
161+
histos.fill(HIST("Lambda/hGenLambda"), particle.pt());
162+
break;
163+
164+
case kLambda0Bar:
165+
histos.fill(HIST("AntiLambda/hGenAntiLambda"), particle.pt());
166+
break;
167+
168+
case kXiMinus:
169+
histos.fill(HIST("Xi/hGenXi"), particle.pt());
170+
break;
171+
172+
case kXiPlusBar:
173+
histos.fill(HIST("AntiXi/hGenAntiXi"), particle.pt());
174+
break;
175+
176+
case kOmegaMinus:
177+
histos.fill(HIST("Omega/hGenOmega"), particle.pt());
178+
break;
179+
180+
case kOmegaPlusBar:
181+
histos.fill(HIST("AntiOmega/hGenAntiOmega"), particle.pt());
182+
break;
183+
184+
default:
185+
break;
186+
}
160187
}
161188

162189
int daughtersIdSlice[2];
@@ -178,13 +205,13 @@ struct OnTheFlyDecayer {
178205
float p = std::sqrt(particle.px() * particle.px() + particle.py() * particle.py() + particle.pz() * particle.pz());
179206
float y; // As https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1943
180207

181-
if ((p - particle.pz()) < VerySmall) {
208+
if ((p - particle.pz()) < Tolerance) {
182209
eta = (particle.pz() < 0.0f) ? -100.0f : 100.0f;
183210
} else {
184211
eta = 0.5f * std::log((p + particle.pz()) / (p - particle.pz()));
185212
}
186213

187-
if ((particle.e() - particle.pz()) < VerySmall) {
214+
if ((particle.e() - particle.pz()) < Tolerance) {
188215
y = (particle.pz() < 0.0f) ? -100.0f : 100.0f;
189216
} else {
190217
y = 0.5f * std::log((particle.e() + particle.pz()) / (particle.e() - particle.pz()));
@@ -196,7 +223,7 @@ struct OnTheFlyDecayer {
196223
particle.flags(), motherSpan, daughtersIdSlice, particle.weight(),
197224
particle.px(), particle.py(), particle.pz(), particle.e(),
198225
particle.vx(), particle.vy(), particle.vz(), particle.vt(),
199-
phi, eta, pt, p, y);
226+
phi, eta, pt, p, y, !canDecay(particle.pdgCode()), true);
200227
}
201228

202229
int daughtersIdSlice[2] = {-1, -1};
@@ -217,26 +244,51 @@ struct OnTheFlyDecayer {
217244
float p = std::sqrt(dau.px() * dau.px() + dau.py() * dau.py() + dau.pz() * dau.pz());
218245
float y; // Conditional as https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1943
219246

220-
if ((p - dau.pz()) < VerySmall) {
247+
if ((p - dau.pz()) < Tolerance) {
221248
eta = (dau.pz() < 0.0f) ? -100.0f : 100.0f;
222249
} else {
223250
eta = 0.5f * std::log((p + dau.pz()) / (p - dau.pz()));
224251
}
225252

226-
if ((dau.e() - dau.pz()) < VerySmall) {
253+
if ((dau.e() - dau.pz()) < Tolerance) {
227254
y = (dau.pz() < 0.0f) ? -100.0f : 100.0f;
228255
} else {
229256
y = 0.5f * std::log((dau.e() + dau.pz()) / (dau.e() - dau.pz()));
230257
}
231258

259+
switch (dau.pdgCode()) {
260+
case kElectron:
261+
histos.fill(HIST("GeneratedElectron/hGenEl"), pt);
262+
break;
263+
264+
case kMuonMinus:
265+
histos.fill(HIST("GeneratedMuon/hGenMu"), pt);
266+
break;
267+
268+
case kPiPlus:
269+
histos.fill(HIST("GeneratedPion/hGenPi"), pt);
270+
break;
271+
272+
case kKPlus:
273+
histos.fill(HIST("GeneratedKaon/hGenKa"), pt);
274+
break;
275+
276+
case kProton:
277+
histos.fill(HIST("GeneratedProton/hGenPr"), pt);
278+
break;
279+
280+
default:
281+
break;
282+
}
283+
232284
// TODO: Particle status code
233285
// TODO: Expression columns
234286
// TODO: vt
235287
tableMcParticlesWithDau(mother.mcCollisionId(), dau.pdgCode(), 1,
236288
mother.flags(), motherSpan, daughtersIdSlice, mother.weight(),
237289
dau.px(), dau.py(), dau.pz(), dau.e(),
238290
dau.vx(), dau.vy(), dau.vz(), mother.vt(),
239-
phi, eta, pt, p, y);
291+
phi, eta, pt, p, y, dau.isAlive(), false);
240292
}
241293
}
242294
}

0 commit comments

Comments
 (0)