|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
|
#include "GeneratorInterface/TauolaInterface/interface/read_particles_from_HepMC.h"
using namespace TauSpinner;
/*******************************************************************************
Get daughters of HepMC::GenParticle
Recursively searches for final-state daughters of 'x'
*******************************************************************************/
inline std::vector<SimpleParticle> *getDaughters(HepMC::GenParticle *x) {
std::vector<SimpleParticle> *daughters = new std::vector<SimpleParticle>();
if (!x->end_vertex())
return daughters;
// Check decay products of 'x'
for (HepMC::GenVertex::particles_out_const_iterator p = x->end_vertex()->particles_out_const_begin();
p != x->end_vertex()->particles_out_const_end();
++p) {
HepMC::GenParticle *pp = *p;
HepMC::FourVector mm = pp->momentum();
// If the daughter of 'x' has its end vertex - recursively read
// all of its daughters.
if (pp->end_vertex() && pp->pdg_id() != 111) {
std::vector<SimpleParticle> *sub_daughters = getDaughters(pp);
daughters->insert(daughters->end(), sub_daughters->begin(), sub_daughters->end());
delete sub_daughters;
}
// Otherwise - add this particle to the list of daughters.
else if (pp->pdg_id() != x->pdg_id()) {
SimpleParticle tp(mm.px(), mm.py(), mm.pz(), mm.e(), pp->pdg_id());
daughters->push_back(tp);
}
}
return daughters;
}
/*******************************************************************************
Find last self
Recursively finds the last particle with the same PDG ID
on the list of its decay products
*******************************************************************************/
inline HepMC::GenParticle *findLastSelf(HepMC::GenParticle *x) {
if (!x->end_vertex())
return x;
for (HepMC::GenVertex::particle_iterator p = x->end_vertex()->particles_begin(HepMC::children);
p != x->end_vertex()->particles_end(HepMC::children);
++p) {
if ((*p)->pdg_id() == x->pdg_id())
return findLastSelf(*p);
}
return x;
}
inline bool isFirst(HepMC::GenParticle *x) {
for (HepMC::GenVertex::particle_iterator p = x->production_vertex()->particles_begin(HepMC::parents);
p != x->production_vertex()->particles_end(HepMC::parents);
++p) {
if (x->pdg_id() == (*p)->pdg_id())
return false;
}
return true;
}
/*******************************************************************************
Read HepMC::GenEvent.
Read HepMC event from data file
and return particles filtered out from the event.
This routine is prepared for use with files generated by Pythia8.
Fills:
'X' - Heavy particle (W+/-, H+/-, H, Z)
'tau' - first tau
'tau2' - second tau or nu_tau, if 'X' decays to one tau only
'tau_daughters' - daughters of 'tau'
'tau2_daughters' - daughters of 'tau2' or empty list, if 'tau2' is nu_tau.
Returns:
0 - no more events to read (finished processing the file)
1 - no decay found in the event (finished processing current event)
2 - decay found and processed correctly.
Event will continue to be processed
with next function call.
*******************************************************************************/
int readParticlesFromHepMC(const HepMC::GenEvent *event,
SimpleParticle &X,
SimpleParticle &tau,
SimpleParticle &tau2,
std::vector<SimpleParticle> &tau_daughters,
std::vector<SimpleParticle> &tau2_daughters) {
if (event == nullptr)
return 1;
// Exctract particles from event
HepMC::GenParticle *hX = nullptr, *hTau = nullptr, *hTau2 = nullptr;
for (HepMC::GenEvent::particle_const_iterator it = event->particles_begin(); it != event->particles_end(); ++it) {
int pdgid = (*it)->pdg_id();
if ((abs(pdgid) == 37 || pdgid == 25 || abs(pdgid) == 24 || pdgid == 23) && (*it)->end_vertex() &&
(*it)->end_vertex()->particles_out_size() >= 2) {
hX = *it;
HepMC::GenParticle *LhX = hX;
if (!isFirst(hX))
continue;
findLastSelf(LhX);
hTau = hTau2 = nullptr;
for (HepMC::GenVertex::particle_iterator it2 = LhX->end_vertex()->particles_begin(HepMC::children);
it2 != LhX->end_vertex()->particles_end(HepMC::children);
++it2) {
if (abs((*it2)->pdg_id()) == 15 && (*it2)->end_vertex()) {
if (!hTau)
hTau = *it2;
else if (!hTau2)
hTau2 = *it2;
else {
std::cout << "TauSpiner: three taus in one decay" << std::endl;
return 1;
}
}
if (abs((*it2)->pdg_id()) == 16) {
if (!hTau2)
hTau2 = *it2;
else {
std::cout << "TauSpiner: two neutrinos or two taus and neutrino in one decay" << std::endl;
return 1;
}
}
}
if (hTau && hTau2)
break;
}
}
if (!hX)
return 1;
if (!hTau || !hTau2) {
//std::cout<<"TauSpiner: boson found but no proper tau pair or tau + neutrino found."<<std::endl;
return 1;
}
// Check for self-decaying taus
hTau = findLastSelf(hTau);
hTau2 = findLastSelf(hTau2);
// Fill SimpleParticles from HepMC particles
X.setPx(hX->momentum().px());
X.setPy(hX->momentum().py());
X.setPz(hX->momentum().pz());
X.setE(hX->momentum().e());
X.setPdgid(hX->pdg_id());
tau.setPx(hTau->momentum().px());
tau.setPy(hTau->momentum().py());
tau.setPz(hTau->momentum().pz());
tau.setE(hTau->momentum().e());
tau.setPdgid(hTau->pdg_id());
tau2.setPx(hTau2->momentum().px());
tau2.setPy(hTau2->momentum().py());
tau2.setPz(hTau2->momentum().pz());
tau2.setE(hTau2->momentum().e());
tau2.setPdgid(hTau2->pdg_id());
// Create list of tau daughters
std::vector<SimpleParticle> *buf = getDaughters(hTau);
tau_daughters.clear();
tau_daughters.insert(tau_daughters.end(), buf->begin(), buf->end());
delete buf;
// Second particle can be 2nd tau. In that case - read its daughters.
// Otherwise it is nu_tau~
buf = getDaughters(hTau2);
tau2_daughters.clear();
tau2_daughters.insert(tau2_daughters.end(), buf->begin(), buf->end());
delete buf;
return 0;
}
|