Execute method.
AV: In case we have 3 particles, we try to add a vertex that corresponds to 1->2 and 1->1 splitting. AV: In case we have 4 particles, we can try to do that as well.
YH: In the case of Sherpa with HEPMC_TREE_LIKE: 1, where the YH: incoming/outgoing particles of the signal process have no YH: production/end vertices, this treatment can produce a loop. YH: Skip it by setting IgnoreSemiDisconnected = True.
Remove loops inside decay chains AV: Please note that this approach would distort some branching ratios. If some particle would have decay products with bad PDG ids, after the operation below the visible branching ratio of these decays would be zero.
Check that all particles coming into the decay vertex of a decay-loop particle candidate came from the same production vertex
36 {
38
41
42 if (!
m_looper.loop_particles().empty() || !
m_looper.loop_vertices().empty()) {
45 ATH_MSG_DEBUG(
"Please use MC::Loops::findLoops for this event to obtain all particles and vertices in the loops");
46 }
47 auto old_momentum =
evt->momentum_unit();
48 auto old_length =
evt->length_unit();
54
56 for (auto ip: *evt) {
57
58 if (!ip) continue;
60 if (newpid ==
m_pidmap.end())
continue;
61 ip->set_pdg_id(newpid->second);
62
65 }
66 }
67
69
70 auto cycles = std::make_shared<HepMC3::IntAttribute>(1);
72 }
73
74
75 if (
evt->weights().empty()) {
76 ATH_MSG_DEBUG(
"Adding a unit weight to empty event weight vector");
77 evt->weights().push_back(1);
78 }
79
82 for (
auto iv:
evt->vertices()) {
83 if (iv->position() == nullpos) {
84 ATH_MSG_DEBUG(
"Setting representative event position vertex");
86 break;
87 }
88 }
89 }
90
91
92 std::vector<HepMC::GenParticlePtr> beams_t;
94 if (
p->status() == 4) beams_t.push_back(p);
95 }
96 if (beams_t.size() > 2) {
97 ATH_MSG_INFO(
"Invalid number of beam particles " << beams_t.size() <<
". Will try to fix.");
98 std::vector<HepMC::GenParticlePtr> bparttoremove;
99 for (const auto& bpart : beams_t) {
100 if (bpart->id() == 0 && bpart->production_vertex()) bparttoremove.push_back(bpart);
101 }
102 for (auto bpart: bparttoremove) {
103 bpart->production_vertex()->remove_particle_out(bpart);
104 }
105 }
106
107
108
109
111 double units_problem = -1.;
112
113 if (beams_t.size() > 0 ){
115 if (
p->momentum().pz() > 1e9) units_problem =
p->momentum().pz();
116 }
117
118 } else {
120 if (p &&
p->momentum().pz() > 1e9) units_problem =
p->momentum().pz();
121 }
122 }
123 if (units_problem>0){
124 ATH_MSG_INFO(
"Apparent units problem; beam particles have z-momentum " << units_problem <<
" in MeV. Will divide by 1000.");
126 }
127 }
128
129
130 std::vector<HepMC::GenParticlePtr> semi_disconnected, decay_loop_particles;
131 for (
auto ip :
evt->particles()) {
132
133 if (!ip) continue;
134 bool particle_to_fix = false;
135 int abspid = std::abs(
ip->pdg_id());
136 auto vProd =
ip->production_vertex();
137 auto vEnd =
ip->end_vertex();
139 if ( (!vProd || vProd->id() == 0) && vEnd &&
ip->status() != 4) {
140 particle_to_fix = true;
141 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without production vertex! HepMC status = " <<
ip->status());
142 }
144 if (vProd && !vEnd &&
ip->status() != 1) {
145 particle_to_fix = true;
146 ATH_MSG_DEBUG(
"Found particle " <<
ip->pdg_id() <<
" without decay vertex! HepMC status = " <<
ip->status());
147 }
148 if (particle_to_fix) semi_disconnected.push_back(ip);
149
150 if (abspid == 43 || abspid == 44 || abspid == 30353 || abspid == 30343) {
151 decay_loop_particles.push_back(std::move(ip));
152 }
153 }
154
158
163 if ( !
m_ignoreSemiDisconnected && (semi_disconnected.size() == 4 || semi_disconnected.size() == 3 || semi_disconnected.size() == 2)) {
164 size_t no_endv = 0;
165 size_t no_prov = 0;
168 for (const auto& part : semi_disconnected) {
169 if (!
part->production_vertex() || !
part->production_vertex()->id()) {
170 no_prov++;
sum +=
part->momentum();
171 }
172 if (!
part->end_vertex()) {
173 no_endv++;
sum -=
part->momentum();
174 }
177 }
178 ATH_MSG_INFO(
"Heuristics: found " << semi_disconnected.size() <<
" semi-disconnected particles. Momentum sum is " << sum <<
" Standalone vertices " <<
standalone.size());
179 bool standalonevertex = (
standalone.size() == 1 && (*
standalone.begin())->particles_in().size() + (*
standalone.begin())->particles_out().size() == semi_disconnected.size());
181 if (! standalonevertex && no_endv && no_prov && ( no_endv + no_prov == semi_disconnected.size() )) {
182 if (std::abs(
sum.px()) < 1e-2 && std::abs(
sum.py()) < 1e-2 && std::abs(
sum.pz()) < 1e-2 ) {
183 ATH_MSG_INFO(
"Try " << no_endv <<
"->" << no_prov <<
" splitting/merging.");
185 for (auto part : semi_disconnected) {
186 if (!
part->production_vertex() ||
part->production_vertex()->id() == 0)
v->add_particle_out(std::move(part));
187 }
188 for (auto part : semi_disconnected) {
189 if (!
part->end_vertex())
v->add_particle_in(std::move(part));
190 }
191 evt->add_vertex(std::move(v));
192 }
193 }
194 }
195
200 for (auto part: decay_loop_particles) {
202 auto vend =
part->end_vertex();
203 auto vprod =
part->production_vertex();
204 if (!vprod || !vend) continue;
205 bool loop_in_decay = true;
208 auto sisters = vend->particles_in();
209 for (auto sister: sisters) {
210 if (vprod != sister->production_vertex()) loop_in_decay = false;
211 }
212 if (!loop_in_decay) continue;
213
215 auto daughters = vend->particles_out();
216 for (auto p : daughters) vprod->add_particle_out(std::move(p));
217 for (auto sister : sisters) {
218 vprod->remove_particle_out(sister);
219 vend->remove_particle_in(sister);
220 evt->remove_particle(std::move(sister));
221 }
222 evt->remove_vertex(std::move(vend));
223
224 }
225
226
227 std::vector<HepMC::GenParticlePtr> toremove;
228 for (
auto ip:
evt->particles()) {
229
230 if (!ip) continue;
232
233 bool bad_particle = false;
234
236 bad_particle = true;
239 if (
msgLvl( MSG::DEBUG ) ) HepMC::Print::line(ip);
240 }
241
243 bad_particle = true;
246 if (
msgLvl( MSG::DEBUG ) )HepMC::Print::line(ip);
247 }
248
249 int abs_pdg_id = std::abs(
ip->pdg_id());
250 bool is_decayed_weak_boson = ( abs_pdg_id == 23 || abs_pdg_id == 24 || abs_pdg_id == 25 ) &&
ip->end_vertex();
252 bad_particle = true;
255 if (
msgLvl( MSG::DEBUG ) ) HepMC::Print::line(ip);
256 }
257
258 if (bad_particle) toremove.push_back(std::move(ip));
259 }
260
261
262 const int num_particles_orig =
evt->particles().size();
263
264
265 if (!toremove.empty()) {
266 ATH_MSG_DEBUG(
"Cleaning event record of " << toremove.size() <<
" bad particles");
267 for (
auto part: toremove)
evt->remove_particle(std::move(part));
268 }
269
271 int purged=0;
272 do {
273 purged=0;
274 const std::vector <HepMC::GenParticlePtr> allParticles=
evt->particles();
275 for(auto p : allParticles) {
277 if(
p->status() == 2 && !end_v) {
278 evt->remove_particle(std::move(p));
279 ++purged;
281 }
282 }
283 }
284 while (purged>0);
285 }
286
287 const int num_particles_filt =
evt->particles().size();
288
289 if(num_particles_orig!=num_particles_filt) {
290
291 ATH_MSG_INFO(
"Particles filtered: " << num_particles_orig <<
" -> " << num_particles_filt);
292 }
293
294 }
295 return StatusCode::SUCCESS;
296}
#define ATH_MSG_WARNING(x)
bool msgLvl(const MSG::Level lvl) const
DataModel_detail::const_iterator< DataVector > const_iterator
bool isPID0(const HepMC::ConstGenParticlePtr &p) const
bool isNonTransportableInDecayChain(const HepMC::ConstGenParticlePtr &p) const
std::map< int, int > m_replacedpid_counts
map of counters of replacements.
MC::Loops< HepMC::GenEvent, HepMC::ConstGenParticlePtr, HepMC::ConstGenVertexPtr > m_looper
member to detect loops
bool isSimpleLoop(const HepMC::ConstGenParticlePtr &p) const
HepMC3::FourVector FourVector
void set_signal_process_vertex(GenEvent *e, T &v)
ConstGenVertexPtr signal_process_vertex(const GenEvent *e)
HepMC3::GenParticlePtr GenParticlePtr
GenVertexPtr newGenVertexPtr(const HepMC3::FourVector &pos=HepMC3::FourVector::ZERO_VECTOR(), const int i=0)
HepMC3::ConstGenVertexPtr ConstGenVertexPtr
HepMC3::GenEvent GenEvent
void MeVToGeV(HepMC::GenEvent *evt)
1.16.1