80{
81
82
83
84
85 static const bool is_sherpa = [this]() {
86 bool is_sherpa = false;
88
89
90
91
93
94 if (
m_metaStore->retrieve(truthMetaData).isSuccess() && !truthMetaData->empty()){
95
96 const std::string gens = truthMetaData->at(0)->generators();
97 is_sherpa = (gens.find("sherpa")==std::string::npos &&
98 gens.find("Sherpa")==std::string::npos &&
99 gens.find("SHERPA")==std::string::npos) ? false : true;
100 }
101 else {
102 ATH_MSG_WARNING(
"Found xAODTruthMetaDataContainer empty! Configuring to be NOT Sherpa.");
103 }
104 ATH_MSG_INFO(
"From metadata configured: Sherpa? " << is_sherpa);
105 } else {
106 ATH_MSG_WARNING(
"Could not find metadata container in storegate; assuming NOT Sherpa");
107 }
108 return is_sherpa;
109 }();
110
111
112 SG::ReadHandle<xAOD::TruthParticleContainer> truthParticles(
m_particlesKey,ctx);
113 if (!truthParticles.isValid()) {
115 return StatusCode::FAILURE;
116 }
117
118
119 SG::WriteHandle<xAOD::TruthParticleContainer> newParticlesWriteHandle(
m_outputParticlesKey, ctx);
120 ATH_CHECK(newParticlesWriteHandle.record(std::make_unique<xAOD::TruthParticleContainer>(),
121 std::make_unique<xAOD::TruthParticleAuxContainer>()));
123
124
125 unsigned int nParticles = truthParticles->size();
127
128
129 SG::ReadDecorHandle<xAOD::TruthParticleContainer, unsigned int > originReadDecor(
m_originReadDecorKey, ctx);
130 SG::ReadDecorHandle<xAOD::TruthParticleContainer, unsigned int > typeReadDecor(
m_typeReadDecorKey, ctx);
131 SG::ReadDecorHandle<xAOD::TruthParticleContainer, unsigned int > outcomeReadDecor(
m_outcomeReadDecorKey, ctx);
133
134
135 SG::WriteDecorHandle<xAOD::TruthParticleContainer, ElementLink<xAOD::TruthParticleContainer> > linkDecorator(
m_linkDecoratorKey, ctx);
136 SG::WriteDecorHandle<xAOD::TruthParticleContainer, unsigned int > originDecorator(
m_originDecoratorKey, ctx);
137 SG::WriteDecorHandle<xAOD::TruthParticleContainer, unsigned int > typeDecorator(
m_typeDecoratorKey, ctx);
138 SG::WriteDecorHandle<xAOD::TruthParticleContainer, unsigned int > outcomeDecorator(
m_outcomeDecoratorKey, ctx);
143
144
146 std::vector<int>
entries = m_parser->evaluateAsVector();
148
149 if (nParticles != nEntries ) {
150 ATH_MSG_ERROR(
"Sizes incompatible! Are you sure your selection string used TruthParticles?");
151 return StatusCode::FAILURE;
152 } else {
153
154
155
156
157
158
159
160 bool SherpaW = false;
161 bool SherpaZ = false;
163 SherpaW = true;
164 }
166 SherpaZ = true;
167 }
168 if ((SherpaW or SherpaZ) && is_sherpa){
170 SherpaW = false;
171 SherpaZ = false;
172 }
173 }
174
175 if ((SherpaW || SherpaZ) && is_sherpa){
176
177 std::vector<const xAOD::TruthParticle*> status_nonPhysical;
178 for (
unsigned int i=0;
i<nParticles; ++
i) {
179
180 if (!truthParticles->at(i)) continue;
181
183
184 if (!
MC::isPhysical(truthParticles->at(i))) status_nonPhysical.push_back( truthParticles->at(i) );
185 }
186
187 if ((status_nonPhysical.size()==2 || status_nonPhysical.size()==4 || status_nonPhysical.size()==6) && (SherpaZ || SherpaW)){
188
189 int gens[3] = {0,0,0};
190 for (
size_t i=0;
i<status_nonPhysical.size();++
i){
191 if (status_nonPhysical[i]->absPdgId()<13) gens[0]++;
192 else if (status_nonPhysical[i]->absPdgId()<15) gens[1]++;
193 else gens[2]++;
194 }
195
196 if (gens[0]>2 || gens[0]==1 || gens[1]>2 || gens[1]==1 || gens[2]>2 || gens[2]==1){
197
198 ATH_MSG_VERBOSE(
"Too many leptons of one generation. Cannot make bosons. Give up");
199 return StatusCode::SUCCESS;
200 }
201 std::vector<const xAOD::TruthParticle*> boson;
202 for (
size_t i=0;
i<status_nonPhysical.size();++
i){
203 if (status_nonPhysical[i]->absPdgId()<13) boson.push_back(status_nonPhysical[i]);
204 else if (gens[0]==0 && status_nonPhysical[i]->absPdgId()<15) boson.push_back(status_nonPhysical[i]);
205 else if (gens[0]==0 && gens[1]==0) boson.push_back(status_nonPhysical[i]);
206 if (boson.size()==2){
207
208 int pdg_id=0;
209
210 if (boson[0]->
pdgId()==-boson[1]->
pdgId()) pdg_id=23;
211 else if (std::abs(boson[0]->
pdgId()+boson[1]->
pdgId())!=1){
212
214 } else {
215
216 pdg_id=24*(boson[0]->pdgId()+boson[1]->
pdgId());
217 }
218 if ( (SherpaW &&
MC::isW(pdg_id)) ||
219 (SherpaZ &&
MC::isZ(pdg_id)) ){
220
222 newParticlesWriteHandle->push_back( xTruthParticle );
224 if (boson[0]->hasProdVtx()) {
225 if ((boson[0]->prodVtx()->nIncomingParticles() > 0) && (boson[0]->prodVtx()->incomingParticle(0)!=nullptr)) {
226 motherIDDecorator(*xTruthParticle) = boson[0]->prodVtx()->incomingParticle(0)->pdgId();
227 } else {motherIDDecorator(*xTruthParticle) = 0;}
228 } else {motherIDDecorator(*xTruthParticle) = 0;}
229 if (boson[0]->hasDecayVtx()) {
230 if ((boson[0]->decayVtx()->nOutgoingParticles() > 0) && (boson[0]->decayVtx()->outgoingParticle(0)!=nullptr)) {
231 daughterIDDecorator(*xTruthParticle) = boson[0]->decayVtx()->outgoingParticle(0)->pdgId();
232 } else {daughterIDDecorator(*xTruthParticle) = 0;}
233 } else {daughterIDDecorator(*xTruthParticle) = 0;}
234 }
235
237
240
242 xTruthParticle->
setM(new_mom.M());
243 xTruthParticle->
setPx(new_mom.Px());
244 xTruthParticle->
setPy(new_mom.Py());
245 xTruthParticle->
setPz(new_mom.Pz());
246 xTruthParticle->
setE(new_mom.E());
247 }
248
249 boson.clear();
250
251 if (gens[0]==0 && gens[1]==0) gens[2]=0;
252 else if (gens[0]==0) gens[1]=0;
253 else gens[0]=0;
254 }
255 }
256 }
257 if (status_nonPhysical.size()==1 || status_nonPhysical.size()==3 || status_nonPhysical.size()==5 || status_nonPhysical.size()>6){
258 ATH_MSG_WARNING(status_nonPhysical.size() <<
" leptons found in the Sherpa event record. Not sure how to deal with this.");
259 }
260 return StatusCode::SUCCESS;
261 }
262 for (
unsigned int i=0;
i<nParticles; ++
i) {
263 ElementLink<xAOD::TruthParticleContainer> eltp(*truthParticles,i);
265
266
269 bool same_as_mother = false;
270 bool same_as_daughter = false;
274 same_as_mother = true;
275 }
276 }
277 }
281 same_as_daughter = true;
282 }
283 }
284 }
285 if (same_as_mother && same_as_daughter){
287 continue;
288 }
289 }
291 newParticlesWriteHandle->push_back( xTruthParticle );
296 } else {motherIDDecorator(*xTruthParticle) = 0;}
297 } else {motherIDDecorator(*xTruthParticle) = 0;}
301 } else {daughterIDDecorator(*xTruthParticle) = 0;}
302 } else {daughterIDDecorator(*xTruthParticle) = 0;}
303 }
304
305 *xTruthParticle=*theParticle;
306
307 typeDecorator(*xTruthParticle) =
308 typeReadDecor(*theParticle);
309
310 originDecorator(*xTruthParticle) =
311 originReadDecor(*theParticle);
312
313 outcomeDecorator(*xTruthParticle) =
314 outcomeReadDecor(*theParticle);
315
316 classificationDecorator(*xTruthParticle) =
317 classificationReadDecor(*theParticle);
318
320 static const SG::ConstAccessor<int> TopHadronOriginFlagAcc("TopHadronOriginFlag");
321 hadronOriginDecorator(*xTruthParticle) =
322 TopHadronOriginFlagAcc.withDefault (*theParticle, 0);
323 }
324
326 }
327 }
328 }
329
331 }
332
333 return StatusCode::SUCCESS;
334}
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_WARNING(x)
ServiceHandle< StoreGateSvc > m_metaStore
Handle on the metadata store for init.
SG::ReadHandleKey< xAOD::TruthParticleContainer > m_particlesKey
Gaudi::Property< bool > m_do_compress
SG::ReadDecorHandleKey< xAOD::TruthParticleContainer > m_originReadDecorKey
std::atomic< unsigned int > m_ntotpart
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_typeDecoratorKey
std::atomic< unsigned int > m_npasspart
SG::ReadDecorHandleKey< xAOD::TruthParticleContainer > m_outcomeReadDecorKey
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_linkDecoratorKey
SG::ReadDecorHandleKey< xAOD::TruthParticleContainer > m_classificationReadDecorKey
Gaudi::Property< std::string > m_partString
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_hadronOriginDecoratorKey
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_originDecoratorKey
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_daughterIDDecoratorKey
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_outcomeDecoratorKey
SG::WriteHandleKey< xAOD::TruthParticleContainer > m_outputParticlesKey
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_classificationDecoratorKey
SG::ReadDecorHandleKey< xAOD::TruthParticleContainer > m_typeReadDecorKey
Gaudi::Property< bool > m_do_sherpa
Gaudi::Property< bool > m_keep_navigation_info
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_motherIDDecoratorKey
TLorentzVector FourMom_t
Definition of the 4-momentum type.
void setPy(float value)
Set the y component of the particle's momentum.
void setUid(int value)
Set unique ID.
void setM(float value)
Also store the mass.
void setE(float value)
Set the energy of the particle.
int pdgId() const
PDG ID code.
const TruthVertex_v1 * decayVtx() const
The decay vertex of this particle.
void setPz(float value)
Set the z component of the particle's momentum.
bool hasProdVtx() const
Check for a production vertex on this particle.
bool hasDecayVtx() const
Check for a decay vertex on this particle.
const TruthVertex_v1 * prodVtx() const
The production vertex of this particle.
void setStatus(int value)
Set status code.
void setPdgId(int pid)
Set PDG ID code.
void setPx(float value)
Set the x component of the particle's momentum.
const TruthParticle_v1 * outgoingParticle(size_t index) const
Get one of the outgoing particles.
const TruthParticle_v1 * incomingParticle(size_t index) const
Get one of the incoming particles.
size_t nOutgoingParticles() const
Get the number of outgoing particles.
size_t nIncomingParticles() const
Get the number of incoming particles.
constexpr int INVALID_PARTICLE_ID
bool isSMLepton(const T &p)
APID: the fourth generation leptons are not standard model leptons.
bool isPhysical(const T &p)
Identify if the particle is physical, i.e. is stable or decayed.
TruthMetaDataContainer_v1 TruthMetaDataContainer
Declare the latest version of the truth vertex container.
TruthParticle_v1 TruthParticle
Typedef to implementation.