-
Notifications
You must be signed in to change notification settings - Fork 494
Expand file tree
/
Copy pathMC.h
More file actions
161 lines (136 loc) · 4.2 KB
/
MC.h
File metadata and controls
161 lines (136 loc) · 4.2 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
// Copyright CERN and copyright holders of ALICE O2. This software is
// distributed under the terms of the GNU General Public License v3 (GPL
// Version 3), copied verbatim in the file "COPYING".
//
// See http://alice-o2.web.cern.ch/license for full licensing information.
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
#ifndef MC_H
#define MC_H
#include "Framework/Logger.h"
#include "TPDGCode.h"
namespace MC
{
bool isStable(int pdg)
{
// Decide whether particle (pdg) is stable
// All ions/nucleons are considered as stable
// Nuclear code is 10LZZZAAAI
if (pdg > 1000000000) {
return true;
}
constexpr int kNstable = 18;
int pdgStable[kNstable] = {
kGamma, // Photon
kElectron, // Electron
kMuonPlus, // Muon
kPiPlus, // Pion
kKPlus, // Kaon
kK0Short, // K0s
kK0Long, // K0l
kProton, // Proton
kNeutron, // Neutron
kLambda0, // Lambda_0
kSigmaMinus, // Sigma Minus
kSigmaPlus, // Sigma Plus
3312, // Xsi Minus
3322, // Xsi
3334, // Omega
kNuE, // Electron Neutrino
kNuMu, // Muon Neutrino
kNuTau // Tau Neutrino
};
for (int i = 0; i < kNstable; i++) {
if (pdg == std::abs(pdgStable[i])) {
return true;
}
}
return false;
}
// Ported from AliRoot AliStack::IsPhysicalPrimary
template <typename TMCParticle, typename TMCParticles>
bool isPhysicalPrimary(TMCParticles& mcParticles, TMCParticle const& particle)
{
// Test if a particle is a physical primary according to the following definition:
// Particles produced in the collision including products of strong and
// electromagnetic decay and excluding feed-down from weak decays of strange
// particles.
LOGF(debug, "isPhysicalPrimary for %d", particle.index());
const int ist = particle.statusCode();
const int pdg = std::abs(particle.pdgCode());
// Initial state particle
// Solution for K0L decayed by Pythia6
// ->
if ((ist > 1) && (pdg != 130) && particle.producedByGenerator()) {
return false;
}
if ((ist > 1) && !particle.producedByGenerator()) {
return false;
}
// <-
if (!isStable(pdg)) {
return false;
}
if (particle.producedByGenerator()) {
// Solution for K0L decayed by Pythia6
// ->
if (particle.mother0() != -1) {
auto mother = mcParticles.iteratorAt(particle.mother0());
if (std::abs(mother.pdgCode()) == 130) {
return false;
}
}
// <-
// check for direct photon in parton shower
// ->
if (pdg == 22 && particle.daughter0() != -1) {
LOGF(debug, "D %d", particle.daughter0());
auto daughter = mcParticles.iteratorAt(particle.daughter0());
if (daughter.pdgCode() == 22) {
return false;
}
}
// <-
return true;
}
// Particle produced during transport
LOGF(debug, "M0 %d %d", particle.producedByGenerator(), particle.mother0());
auto mother = mcParticles.iteratorAt(particle.mother0());
int mpdg = std::abs(mother.pdgCode());
// Check for Sigma0
if ((mpdg == 3212) && mother.producedByGenerator()) {
return true;
}
// Check if it comes from a pi0 decay
if ((mpdg == kPi0) && mother.producedByGenerator()) {
return true;
}
// Check if this is a heavy flavor decay product
int mfl = int(mpdg / std::pow(10, int(std::log10(mpdg))));
// Light hadron
if (mfl < 4) {
return false;
}
// Heavy flavor hadron produced by generator
if (mother.producedByGenerator()) {
return true;
}
// To be sure that heavy flavor has not been produced in a secondary interaction
// Loop back to the generated mother
LOGF(debug, "M0 %d %d", mother.producedByGenerator(), mother.mother0());
while (mother.mother0() != -1 && !mother.producedByGenerator()) {
mother = mcParticles.iteratorAt(mother.mother0());
LOGF(debug, "M+ %d %d", mother.producedByGenerator(), mother.mother0());
mpdg = std::abs(mother.pdgCode());
mfl = int(mpdg / std::pow(10, int(std::log10(mpdg))));
}
if (mfl < 4) {
return false;
} else {
return true;
}
}
}; // namespace MC
#endif