75{
76
77
78
79
80 static const bool is_sherpa = [this]() {
81 bool is_sherpa = false;
83
84
85
86
88
89 if (
m_metaStore->retrieve(truthMetaData).isSuccess() && !truthMetaData->empty()){
90
91 const std::string gens = truthMetaData->at(0)->generators();
92 is_sherpa = (gens.find("sherpa")==std::string::npos &&
93 gens.find("Sherpa")==std::string::npos &&
94 gens.find("SHERPA")==std::string::npos) ? false : true;
95 }
96 else {
97 ATH_MSG_WARNING(
"Found xAODTruthMetaDataContainer empty! Configuring to be NOT Sherpa.");
98 }
99 ATH_MSG_INFO(
"From metadata configured: Sherpa? " << is_sherpa);
100 } else {
101 ATH_MSG_WARNING(
"Could not find metadata container in storegate; assuming NOT Sherpa");
102 }
103 return is_sherpa;
104 }();
105
106
107 SG::ReadHandle<xAOD::TruthParticleContainer> truthParticles(
m_particlesKey,ctx);
108 if (!truthParticles.isValid()) {
110 return StatusCode::FAILURE;
111 }
112
113
114 SG::WriteHandle<xAOD::TruthParticleContainer> newParticlesWriteHandle(
m_outputParticlesKey, ctx);
115 ATH_CHECK(newParticlesWriteHandle.record(std::make_unique<xAOD::TruthParticleContainer>(),
116 std::make_unique<xAOD::TruthParticleAuxContainer>()));
118
119
120 unsigned int nParticles = truthParticles->size();
122
123
124 SG::ReadDecorHandle<xAOD::TruthParticleContainer, unsigned int > originReadDecor(
m_originReadDecorKey, ctx);
125 SG::ReadDecorHandle<xAOD::TruthParticleContainer, unsigned int > typeReadDecor(
m_typeReadDecorKey, ctx);
126 SG::ReadDecorHandle<xAOD::TruthParticleContainer, unsigned int > outcomeReadDecor(
m_outcomeReadDecorKey, ctx);
128
129
130 SG::WriteDecorHandle<xAOD::TruthParticleContainer, ElementLink<xAOD::TruthParticleContainer> > linkDecorator(
m_linkDecoratorKey, ctx);
131 SG::WriteDecorHandle<xAOD::TruthParticleContainer, unsigned int > originDecorator(
m_originDecoratorKey, ctx);
132 SG::WriteDecorHandle<xAOD::TruthParticleContainer, unsigned int > typeDecorator(
m_typeDecoratorKey, ctx);
133 SG::WriteDecorHandle<xAOD::TruthParticleContainer, unsigned int > outcomeDecorator(
m_outcomeDecoratorKey, ctx);
138
139
142
143 if (nParticles != nEntries ) {
144 ATH_MSG_ERROR(
"Sizes incompatible! Are you sure your selection string used TruthParticles?");
145 return StatusCode::FAILURE;
146 } else {
147
148
149
150
151
152
153
154 bool SherpaW = false;
155 bool SherpaZ = false;
157 SherpaW = true;
158 }
160 SherpaZ = true;
161 }
162 if ((SherpaW or SherpaZ) && is_sherpa){
164 SherpaW = false;
165 SherpaZ = false;
166 }
167 }
168
169 if ((SherpaW || SherpaZ) && is_sherpa){
170
171 std::vector<const xAOD::TruthParticle*> status_nonPhysical;
172 for (
unsigned int i=0;
i<nParticles; ++
i) {
173
174 if (!truthParticles->at(i)) continue;
175
177
178 if (!
MC::isPhysical(truthParticles->at(i))) status_nonPhysical.push_back( truthParticles->at(i) );
179 }
180
181 if ((status_nonPhysical.size()==2 || status_nonPhysical.size()==4 || status_nonPhysical.size()==6) && (SherpaZ || SherpaW)){
182
183 int gens[3] = {0,0,0};
184 for (
size_t i=0;
i<status_nonPhysical.size();++
i){
185 if (status_nonPhysical[i]->absPdgId()<13) gens[0]++;
186 else if (status_nonPhysical[i]->absPdgId()<15) gens[1]++;
187 else gens[2]++;
188 }
189
190 if (gens[0]>2 || gens[0]==1 || gens[1]>2 || gens[1]==1 || gens[2]>2 || gens[2]==1){
191
192 ATH_MSG_VERBOSE(
"Too many leptons of one generation. Cannot make bosons. Give up");
193 return StatusCode::SUCCESS;
194 }
195 std::vector<const xAOD::TruthParticle*> boson;
196 for (
size_t i=0;
i<status_nonPhysical.size();++
i){
197 if (status_nonPhysical[i]->absPdgId()<13) boson.push_back(status_nonPhysical[i]);
198 else if (gens[0]==0 && status_nonPhysical[i]->absPdgId()<15) boson.push_back(status_nonPhysical[i]);
199 else if (gens[0]==0 && gens[1]==0) boson.push_back(status_nonPhysical[i]);
200 if (boson.size()==2){
201
202 int pdg_id=0;
203
204 if (boson[0]->
pdgId()==-boson[1]->
pdgId()) pdg_id=23;
205 else if (std::abs(boson[0]->
pdgId()+boson[1]->
pdgId())!=1){
206
208 } else {
209
210 pdg_id=24*(boson[0]->pdgId()+boson[1]->
pdgId());
211 }
212 if ( (SherpaW &&
MC::isW(pdg_id)) ||
213 (SherpaZ &&
MC::isZ(pdg_id)) ){
214
216 newParticlesWriteHandle->push_back( xTruthParticle );
218 if (boson[0]->hasProdVtx()) {
219 if ((boson[0]->prodVtx()->nIncomingParticles() > 0) && (boson[0]->prodVtx()->incomingParticle(0)!=nullptr)) {
220 motherIDDecorator(*xTruthParticle) = boson[0]->prodVtx()->incomingParticle(0)->pdgId();
221 } else {motherIDDecorator(*xTruthParticle) = 0;}
222 } else {motherIDDecorator(*xTruthParticle) = 0;}
223 if (boson[0]->hasDecayVtx()) {
224 if ((boson[0]->decayVtx()->nOutgoingParticles() > 0) && (boson[0]->decayVtx()->outgoingParticle(0)!=nullptr)) {
225 daughterIDDecorator(*xTruthParticle) = boson[0]->decayVtx()->outgoingParticle(0)->pdgId();
226 } else {daughterIDDecorator(*xTruthParticle) = 0;}
227 } else {daughterIDDecorator(*xTruthParticle) = 0;}
228 }
229
231
234
236 xTruthParticle->
setM(new_mom.M());
237 xTruthParticle->
setPx(new_mom.Px());
238 xTruthParticle->
setPy(new_mom.Py());
239 xTruthParticle->
setPz(new_mom.Pz());
240 xTruthParticle->
setE(new_mom.E());
241 }
242
243 boson.clear();
244
245 if (gens[0]==0 && gens[1]==0) gens[2]=0;
246 else if (gens[0]==0) gens[1]=0;
247 else gens[0]=0;
248 }
249 }
250 }
251 if (status_nonPhysical.size()==1 || status_nonPhysical.size()==3 || status_nonPhysical.size()==5 || status_nonPhysical.size()>6){
252 ATH_MSG_WARNING(status_nonPhysical.size() <<
" leptons found in the Sherpa event record. Not sure how to deal with this.");
253 }
254 return StatusCode::SUCCESS;
255 }
256 for (
unsigned int i=0;
i<nParticles; ++
i) {
257 ElementLink<xAOD::TruthParticleContainer> eltp(*truthParticles,i);
259
260
263 bool same_as_mother = false;
264 bool same_as_daughter = false;
268 same_as_mother = true;
269 }
270 }
271 }
275 same_as_daughter = true;
276 }
277 }
278 }
279 if (same_as_mother && same_as_daughter){
281 continue;
282 }
283 }
285 newParticlesWriteHandle->push_back( xTruthParticle );
290 } else {motherIDDecorator(*xTruthParticle) = 0;}
291 } else {motherIDDecorator(*xTruthParticle) = 0;}
295 } else {daughterIDDecorator(*xTruthParticle) = 0;}
296 } else {daughterIDDecorator(*xTruthParticle) = 0;}
297 }
298
299 *xTruthParticle=*theParticle;
300
301 typeDecorator(*xTruthParticle) =
302 typeReadDecor(*theParticle);
303
304 originDecorator(*xTruthParticle) =
305 originReadDecor(*theParticle);
306
307 outcomeDecorator(*xTruthParticle) =
308 outcomeReadDecor(*theParticle);
309
310 classificationDecorator(*xTruthParticle) =
311 classificationReadDecor(*theParticle);
312
314 static const SG::ConstAccessor<int> TopHadronOriginFlagAcc("TopHadronOriginFlag");
315 hadronOriginDecorator(*xTruthParticle) =
316 TopHadronOriginFlagAcc.withDefault (*theParticle, 0);
317 }
318
320 }
321 }
322 }
323
325
326 return StatusCode::SUCCESS;
327}
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_WARNING(x)
SG::WriteHandleKey< xAOD::TruthParticleContainer > m_outputParticlesKey
std::atomic< unsigned int > m_ntotpart
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_hadronOriginDecoratorKey
virtual std::vector< int > updateMask(const xAOD::TruthParticleContainer *) const
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_linkDecoratorKey
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_outcomeDecoratorKey
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_classificationDecoratorKey
std::atomic< unsigned int > m_npasspart
SG::ReadDecorHandleKey< xAOD::TruthParticleContainer > m_classificationReadDecorKey
Gaudi::Property< bool > m_buildingW
Gaudi::Property< bool > m_do_sherpa
Gaudi::Property< bool > m_do_compress
SG::ReadDecorHandleKey< xAOD::TruthParticleContainer > m_outcomeReadDecorKey
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_motherIDDecoratorKey
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_daughterIDDecoratorKey
Gaudi::Property< bool > m_keep_navigation_info
Gaudi::Property< bool > m_buildingZ
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_originDecoratorKey
SG::ReadDecorHandleKey< xAOD::TruthParticleContainer > m_typeReadDecorKey
ServiceHandle< StoreGateSvc > m_metaStore
Handle on the metadata store for init.
SG::ReadDecorHandleKey< xAOD::TruthParticleContainer > m_originReadDecorKey
SG::ReadHandleKey< xAOD::TruthParticleContainer > m_particlesKey
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_typeDecoratorKey
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.