18 {
20 int topbar = 0;
21 bool isLastTop = false;
22 bool isFirstTop = false;
23 double topPairInvariantMass = 0.;
24 std::vector<HepMC::ConstGenParticlePtr> tops;
25 std::vector<HepMC::ConstGenVertexPtr> top_vtxs;
26
28 const HepMC::GenEvent* genEvt = (*itr);
29#ifdef HEPMC3
30 for (const auto& mcpart: *genEvt) {
31
32 isLastTop = false;
33
34 isFirstTop = false;
35
36
37 if (std::abs(mcpart->pdg_id()) == 6) {
38
39 isLastTop=true;
40 auto decayVtx = mcpart->end_vertex();
41
42 if (!decayVtx) {
43 ATH_MSG_WARNING(
"top particle" << mcpart <<
" has no valid decay vertex. ");
44 ATH_MSG_WARNING(
"It looks like a Pythia history particle. Skip this particle ");
45 continue;
46 }
47
48
50 for (const auto& child_mcpart: decayVtx->particles_out() ) {
51 if (std::abs(child_mcpart->pdg_id()) == 6) {
52
53 isLastTop = false;
54 break;
55 }
56 }
57
58
59 if (isLastTop) {
60 ATH_MSG_DEBUG(
"Top particle " << mcpart <<
" is found and stored");
61 tops.push_back(mcpart);
62 }
63
64
65
66 if (isLastTop) {
67
68 isFirstTop = false;
69
70 auto prodVtx = mcpart->production_vertex();
71 if (!prodVtx) {
72 ATH_MSG_WARNING(
"Top particle " << mcpart <<
" has no valid production vertex");
73
74 top_vtxs.emplace_back(nullptr);
75 } else {
76
77 while (!isFirstTop && prodVtx) {
78
79 for (const auto& mother_mcpart: prodVtx->particles_in()) {
80
81 if (mother_mcpart->pdg_id() == 6) {
82 isFirstTop = false;
83 prodVtx = mother_mcpart->production_vertex();
84 if (!prodVtx) {
85 ATH_MSG_WARNING(
"mother particle is still a top " << mcpart <<
", but has no valid production vertex");
86 top_vtxs.emplace_back(nullptr);
87 }
88 break;
89 } else {
90
91
92
93 isFirstTop = true;
94 }
95 }
96 }
97
98
99
100 if (isFirstTop) {
101 top_vtxs.push_back(std::move(prodVtx));
102 }
103 }
104 }
105
106
107
108 if (isLastTop) {
109
110 if (tops[tops.size()-1]->pdg_id() == 6)
top++;
else topbar++;
111
112
113
114 if (
top==1 && topbar==1)
break;
115
116
117
118
119
120
121 if (
top > 1 || topbar > 1)
ATH_MSG_INFO(
"More than one top-pair exist in the event.");
122
123
124
125
126 if (
top==2 && topbar==2)
break;
127 }
128 }
129 }
130#else
131 HepMC::GenEvent::particle_const_iterator pitr = genEvt->particles_begin();
132 for (; pitr!=genEvt->particles_end(); ++pitr) {
133
134 isLastTop = false;
135
136 isFirstTop = false;
137 HepMC::GenParticle* mcpart = (*pitr);
138
139
140 if (std::abs(mcpart->pdg_id()) == 6) {
141
142 isLastTop=true;
143 HepMC::GenVertex* decayVtx = mcpart->end_vertex();
144
145 if (!decayVtx) {
146 ATH_MSG_WARNING(
"top particle " << mcpart <<
" has no valid decay vertex. ");
147 ATH_MSG_WARNING(
"It looks like a Pythia history particle. Skip this particle ");
148 continue;
149 }
150
151
153
154 HepMC::GenVertex::particles_out_const_iterator child_mcpartItr = decayVtx->particles_out_const_begin();
155 HepMC::GenVertex::particles_out_const_iterator child_mcpartItrE = decayVtx->particles_out_const_end();
156 for (; child_mcpartItr!=child_mcpartItrE; ++child_mcpartItr) {
157 HepMC::GenParticle* child_mcpart = (*child_mcpartItr);
158 if (std::abs(child_mcpart->pdg_id()) == 6) {
159
160 isLastTop = false;
161 break;
162 }
163 }
164
165
166 if (isLastTop) {
167 ATH_MSG_DEBUG(
"Top particle " << mcpart <<
" is found and stored");
168 tops.push_back(mcpart);
169 }
170
171
172
173 if (isLastTop) {
174
175 isFirstTop = false;
176
177 HepMC::GenVertex* prodVtx = mcpart->production_vertex();
178 if (!prodVtx) {
179 ATH_MSG_WARNING(
"Top particle " << mcpart <<
" has no valid production vertex");
180
181 top_vtxs.push_back(NULL);
182 } else {
183
184 while (!isFirstTop && prodVtx) {
185
186 HepMC::GenVertex::particles_in_const_iterator mother_mcpartItr = prodVtx->particles_in_const_begin();
187 HepMC::GenVertex::particles_in_const_iterator mother_mcpartItrE = prodVtx->particles_in_const_end();
188 for (; mother_mcpartItr != mother_mcpartItrE; ++mother_mcpartItr) {
189 HepMC::GenParticle* mother_mcpart = (*mother_mcpartItr);
190
191 if (mother_mcpart->pdg_id() == 6) {
192 isFirstTop = false;
193 prodVtx = mother_mcpart->production_vertex();
194 if (!prodVtx) {
195 ATH_MSG_WARNING(
"mother particle is still a top " << mcpart <<
", but has no valid production vertex");
196 top_vtxs.push_back(NULL);
197 }
198 break;
199 } else {
200
201
202
203 isFirstTop = true;
204 }
205 }
206 }
207
208
209
210 if (isFirstTop) {
211 top_vtxs.push_back(prodVtx);
212 }
213 }
214 }
215
216
217
218 if (isLastTop) {
219
220 if (tops[tops.size()-1]->pdg_id() == 6)
top++;
else topbar++;
221
222
223
224 if (
top==1 && topbar==1)
break;
225
226
227
228
229
230
231 if (
top > 1 || topbar > 1)
ATH_MSG_INFO(
"More than one top-pair exist in the event.");
232
233
234
235
236 if (
top==2 && topbar==2)
break;
237 }
238
239 }
240 }
241#endif
242 }
243
244
245 if (tops.size()==2) {
246 HepMC::FourVector topSumMomentum(0.,0.,0.,0.);
251 topPairInvariantMass = topSumMomentum.m();
252 }
253
254
255
256
257
258 else if (tops.size() == 4 && top_vtxs.size() == 4) {
259 HepMC::FourVector topSumMomentum_1(0.,0.,0.,0.);
260 HepMC::FourVector topSumMomentum_2(0.,0.,0.,0.);
261 double topPairInvMass_1=0.;
262 double topPairInvMass_2=0.;
263
264
265 auto prodVtx = top_vtxs[0];
266 int top_11 = 0;
267 int top_12 = -1;
268 int top_21 = -1;
269 int top_22 = -1;
270 for (
size_t i = 1;
i < top_vtxs.size(); ++
i) {
271 if (top_vtxs[i] == prodVtx) top_12=
i;
272 else {
273 if (top_21 < 0) top_21 =
i;
275 }
276 }
277
278
279 if (top_12 < 0) {
280 ATH_MSG_ERROR(
"First top particle doesn't share the production vertex to any other top particles. Event failed the filter");
281 setFilterPassed(false);
282 return StatusCode::SUCCESS;
283 }
284 if ((top_21 < 0) or (top_22 < 0)) {
286 setFilterPassed(false);
287 return StatusCode::SUCCESS;
288 }
289
290
291 if (top_vtxs[top_21] != top_vtxs[top_22]) {
292 ATH_MSG_ERROR(
"Production vertex for the second top-pair particles is not the same. Event failed the filter");
293 setFilterPassed(false);
294 return StatusCode::SUCCESS;
295 }
296
297
298 if (tops[top_11]->pdg_id() == tops[top_12]->pdg_id())
ATH_MSG_ERROR(
"first top-pair particles have the same charge");
299 else {
300
305 topPairInvMass_1 = topSumMomentum_1.m();
306 }
307
308
309 if (tops[top_21]->pdg_id() == tops[top_21]->pdg_id())
ATH_MSG_ERROR(
"second top-pair particles have the same charge");
310 else {
311
316 topPairInvMass_2 = topSumMomentum_2.m();
317 }
318
319
320 topPairInvariantMass = std::max(topPairInvMass_1, topPairInvMass_2);
321 }
322
323
324 else ATH_MSG_ERROR(
"Size of the collected final top particles vector " << tops.size() <<
325 " doesn't match to the size of there production vertices vector size " << top_vtxs.size() <<
326 ". Event fails the filter");
327
328
329
330
331
332 if (tops.size() != 2 && tops.size() != 4) {
333 ATH_MSG_WARNING(
"1 or 2 top-pairs are expected in the event. But there are " << tops.size() <<
" final-status top particles.");
334 }
335
336
338 ATH_MSG_DEBUG(
"The top-pair invariant mass is " << topPairInvariantMass <<
" CLHEP::MeV");
339 return StatusCode::SUCCESS;
340}
DataModel_detail::const_iterator< DataVector > const_iterator