18 {
19 int N_quark_t = 0;
20 int N_quark_tbar = 0;
21 int N_quark_t_all = 0;
22 int N_quark_tbar_all = 0;
23 int N_pt_above_cut = 0;
24 int N_pt_above_cut_plus = 0;
25 int N_pt_above_cut_minus = 0;
26
27 int foundlepton[6] = {0, 0, 0, 0, 0, 0};
28 int count_found_leptons = 1;
29 if(m_fourTopsFilter) count_found_leptons = 2;
30
32 const HepMC::GenEvent* genEvt = (*itr);
33#ifdef HEPMC3
34 for (const auto& pitr: *genEvt) {
35 if (std::abs(pitr->pdg_id()) != 6) continue;
36 if ( pitr->pdg_id() == 6 ) N_quark_t_all++;
37 if ( pitr->pdg_id() == -6 ) N_quark_tbar_all++;
38 auto decayVtx = pitr->end_vertex();
39
40 if (!decayVtx) continue;
41
42 if (decayVtx->particles_out().size() < 2) continue;
43 for (const auto& child_mcpart: decayVtx->particles_out()) {
44
45 if (std::abs(child_mcpart->pdg_id()) != 24) continue;
46 if ( pitr->pdg_id() == 6 ) N_quark_t++;
47 if ( pitr->pdg_id() == -6 ) N_quark_tbar++;
48
49 bool useNextVertex = false;
50 auto w_decayVtx = child_mcpart->end_vertex();
51 while (w_decayVtx) {
52 useNextVertex = false;
53 for (const auto& grandchild_mcpart: *w_decayVtx) {
54 int grandchild_pid = grandchild_mcpart->pdg_id();
55 ATH_MSG_DEBUG(
"W (t/tbar) has " << w_decayVtx->particles_out().size() <<
" children and the pdg_id of the next is " << grandchild_pid);
56
57 if (std::abs(grandchild_pid) == 24) {
58 w_decayVtx = grandchild_mcpart->end_vertex();
59
60 if (!w_decayVtx) {
62 break;
63 }
64 useNextVertex = true;
65 break;
66 }
67
68
69 if (grandchild_pid == -11 && foundlepton[0] < count_found_leptons) {
70 if (grandchild_mcpart->momentum().perp() >= m_Ptmin) {
71 foundlepton[0]++;
72 N_pt_above_cut++;
73 N_pt_above_cut_minus++;
74 }
75 }
76 if (grandchild_pid == 11 && foundlepton[1] < count_found_leptons) {
77 if (grandchild_mcpart->momentum().perp() >= m_Ptmin) {
78 foundlepton[1]++;
79 N_pt_above_cut++;
80 N_pt_above_cut_plus++;
81
82 }
83 }
84 if (grandchild_pid == -13 && foundlepton[2] < count_found_leptons) {
85 if (grandchild_mcpart->momentum().perp() >= m_Ptmin) {
86 foundlepton[2]++;
87 N_pt_above_cut++;
88 N_pt_above_cut_minus++;
89 }
90 }
91 if (grandchild_pid == 13 && foundlepton[3] < count_found_leptons) {
92 if (grandchild_mcpart->momentum().perp() >= m_Ptmin) {
93 foundlepton[3]++;
94 N_pt_above_cut++;
95 N_pt_above_cut_plus++;
96 }
97 }
98 if (grandchild_pid == -15 && foundlepton[4] < count_found_leptons) {
99 if (grandchild_mcpart->momentum().perp() >= m_Ptmin) {
100 foundlepton[4]++;
101 N_pt_above_cut++;
102 N_pt_above_cut_minus++;
103 }
104 }
105 if (grandchild_pid == 15 && foundlepton[5] < count_found_leptons) {
106 if (grandchild_mcpart->momentum().perp() >= m_Ptmin) {
107 foundlepton[5]++;
108 N_pt_above_cut++;
109 N_pt_above_cut_plus++;
110 }
111 }
112 }
113
114 if (!useNextVertex) break;
115 }
116 }
117 }
118#else
119 for (HepMC::GenEvent::particle_const_iterator pitr = genEvt->particles_begin(); pitr != genEvt->particles_end(); ++pitr) {
120 if (std::abs((*pitr)->pdg_id()) == 6) {
121 if ( (*pitr)->pdg_id() == 6 ) N_quark_t_all++;
122 if ( (*pitr)->pdg_id() == -6 ) N_quark_tbar_all++;
123
124 int n_daughters = 0;
125
126 HepMC::GenParticle * mcpart = (*pitr);
127 const HepMC::GenVertex * decayVtx = mcpart->end_vertex();
128
129
130 if (decayVtx != 0) n_daughters = decayVtx->particles_out_size();
131
132
133 if (n_daughters >= 2) {
134 HepMC::GenVertex::particles_in_const_iterator child_mcpartItr = decayVtx->particles_out_const_begin();
135 HepMC::GenVertex::particles_in_const_iterator child_mcpartItrE = decayVtx->particles_out_const_end();
136 for (; child_mcpartItr != child_mcpartItrE; ++child_mcpartItr) {
137 HepMC::GenParticle * child_mcpart = (*child_mcpartItr);
138
139
140 if (std::abs(child_mcpart->pdg_id()) == 24) {
141 if ( (*pitr)->pdg_id() == 6 ) N_quark_t++;
142 if ( (*pitr)->pdg_id() == -6 ) N_quark_tbar++;
143
144 bool useNextVertex = false;
145 const HepMC::GenVertex * w_decayVtx = child_mcpart->end_vertex();
146
147 while (w_decayVtx) {
148
149 useNextVertex = false;
150 int mcpart_n_particles_out = w_decayVtx->particles_out_size();
151
152 HepMC::GenVertex::particles_out_const_iterator grandchild_mcpartItr = w_decayVtx->particles_out_const_begin();
153 HepMC::GenVertex::particles_out_const_iterator grandchild_mcpartItrE = w_decayVtx->particles_out_const_end();
154
155 for (; grandchild_mcpartItr != grandchild_mcpartItrE; ++grandchild_mcpartItr) {
156 HepMC::GenParticle * grandchild_mcpart = (*grandchild_mcpartItr);
157 int grandchild_pid = grandchild_mcpart->pdg_id();
158
159 ATH_MSG_DEBUG(
"W (t/tbar) has " << mcpart_n_particles_out <<
" children and the pdg_id of the next is " << grandchild_pid);
160
161
162 if (std::abs(grandchild_pid) == 24) {
163 w_decayVtx = grandchild_mcpart->end_vertex();
164
165
166 if (!w_decayVtx) {
168 break;
169 }
170
171 useNextVertex = true;
172 break;
173 }
174
175
176 if (grandchild_pid == -11 && foundlepton[0] < count_found_leptons) {
177 if (grandchild_mcpart->momentum().perp() >= m_Ptmin) {
178 foundlepton[0]++;
179 N_pt_above_cut++;
180 N_pt_above_cut_minus++;
181 }
182 }
183 if (grandchild_pid == 11 && foundlepton[1] < count_found_leptons) {
184 if (grandchild_mcpart->momentum().perp() >= m_Ptmin) {
185 foundlepton[1]++;
186 N_pt_above_cut++;
187 N_pt_above_cut_plus++;
188
189 }
190 }
191 if (grandchild_pid == -13 && foundlepton[2] < count_found_leptons) {
192 if (grandchild_mcpart->momentum().perp() >= m_Ptmin) {
193 foundlepton[2]++;
194 N_pt_above_cut++;
195 N_pt_above_cut_minus++;
196 }
197 }
198 if (grandchild_pid == 13 && foundlepton[3] < count_found_leptons) {
199 if (grandchild_mcpart->momentum().perp() >= m_Ptmin) {
200 foundlepton[3]++;
201 N_pt_above_cut++;
202 N_pt_above_cut_plus++;
203 }
204 }
205 if (grandchild_pid == -15 && foundlepton[4] < count_found_leptons) {
206 if (grandchild_mcpart->momentum().perp() >= m_Ptmin) {
207 foundlepton[4]++;
208 N_pt_above_cut++;
209 N_pt_above_cut_minus++;
210 }
211 }
212 if (grandchild_pid == 15 && foundlepton[5] < count_found_leptons) {
213 if (grandchild_mcpart->momentum().perp() >= m_Ptmin) {
214 foundlepton[5]++;
215 N_pt_above_cut++;
216 N_pt_above_cut_plus++;
217 }
218 }
219
220 }
221
222
223 if (!useNextVertex) break;
224 }
225 }
226 }
227 }
228 }
229 }
230#endif
231 }
232
233 ATH_MSG_INFO(
"Found " << N_quark_t_all <<
" t quarks in event record");
234 ATH_MSG_INFO(
"Found " << N_quark_tbar_all <<
" tbar quarks in event record");
235 ATH_MSG_INFO(
"Found " << N_quark_t <<
" t -> W X decays");
236 ATH_MSG_INFO(
"Found " << N_quark_tbar <<
" tbar -> W X decays");
237 ATH_MSG_INFO(
"Num leptons from W decays from tops with lepton pt above " << m_Ptmin/1000.0 <<
" Gaudi::Units::GeV/c = " << N_pt_above_cut);
238
239 int count_tops = 1;
240 if(m_fourTopsFilter) count_tops = 2;
241 if (N_quark_t_all < count_tops || N_quark_tbar_all < count_tops) {
242
243 ATH_MSG_ERROR(
"No t or tbar quarks were found in a (presumably) ttbar event! Event is rejected.");
244 setFilterPassed(false);
245 return StatusCode::SUCCESS;
246 }
247
248 if (N_quark_t < count_tops || N_quark_tbar < count_tops) {
249
250 ATH_MSG_ERROR(
"No t or tbar quarks were found decaying to W in a (presumably) ttbar event! Event is rejected. Event dump follows.");
251 int event = 0;
253 event++;
254 const HepMC::GenEvent* genEvt = (*itr);
256 for (const auto& mcpart: *genEvt) {
258 int pid = mcpart->pdg_id();
259 ATH_MSG_ERROR(
"In event (from MC collection) " << event <<
" particle number " << part <<
" has pdg_id = " << pid);
260
261 auto decayVtx = mcpart->end_vertex();
262
263 if ( !decayVtx ) continue;
264 int part_child=0;
265 for (const auto& child_mcpart: *decayVtx) {
266 part_child++;
267 int child_pid = child_mcpart->pdg_id();
268 ATH_MSG_ERROR(
" child " << part_child <<
" with pdg_id = " << child_pid);
269 }
270 }
271 }
272 setFilterPassed(false);
273 return StatusCode::SUCCESS;
274 }
275
276
277 if (m_numLeptons < 0) {
278 if ( (N_quark_t >= 2 || N_quark_tbar >= 2) && !m_fourTopsFilter) {
279 ATH_MSG_WARNING(
"More than one t -> W X or tbar -> W X decays found. Event is accepted anyway.");
280 }
281 if ( (N_quark_t > 2 || N_quark_tbar > 2) && m_fourTopsFilter) {
282 ATH_MSG_WARNING(
"More than one t -> W X or tbar -> W X decays found. Event is accepted anyway.");
283 }
284 setFilterPassed(N_pt_above_cut > 0);
285 } else {
286 if(m_fourTopsFilter){
287 if(m_SSMLFilter) setFilterPassed( (N_pt_above_cut >= m_numLeptons) && (N_pt_above_cut_plus >= 2 || N_pt_above_cut_minus >= 2));
288 else setFilterPassed(N_pt_above_cut >= m_numLeptons);}
289 else setFilterPassed(N_pt_above_cut == m_numLeptons);
290 }
291
292 return StatusCode::SUCCESS;
293}
DataModel_detail::const_iterator< DataVector > const_iterator