62 {
63
65 const std::vector<Trk::VxTrackAtVertex *> * tracks =
vertex->vxTrackAtVertex();
66 if ( tracks !=nullptr ) {
67
68
71 if (
sc.isFailure() ) {
74 msg() <<
"Zero pointer returned." <<
endmsg;
75 }
76 return nullptr;
77 }
78
79
81 const HepMC::GenEvent* genEvent= ( *it );
82
83 if( genEvent->vertices_empty() ) {
85 return nullptr;
86 }
87#ifdef HEPMC3
88 auto pv = genEvent->vertices()[0];
89#else
90 auto pv = *(genEvent->vertices_begin());
91#endif
92
93
94
95
96
97 const auto& pv_pos =
pv ->position();
98 double pv_r = pv_pos.perp();
99 double pv_z = pv_pos.z();
100
101
102
103 std::map<int,HepMC::ConstGenVertexPtr> vertex_ids;
104#ifdef HEPMC3
105 for (const auto& vtx: genEvent->vertices()){
106 const auto& lv_pos = vtx->position();
107 if ( std::abs ( lv_pos.perp() - pv_r ) <
m_r_tol && std::abs ( lv_pos.z() - pv_z ) <
m_z_tol ) {vertex_ids[vtx->id()] = vtx;}
108 }
109#else
110 for ( HepMC::GenEvent::vertex_const_iterator i = genEvent->vertices_begin(); i != genEvent->vertices_end() ;++i ) {
112 const auto& lv_pos = vtx->position();
114 }
115#endif
116
117
118
121
122 if (
sc.isFailure() ) {
125 msg() <<
"Zero pointer returned" <<
endmsg;
126 }
127 return nullptr;
128 }
129
130
131 std::vector<Trk::VxTrackAtVertex *>::const_iterator
vt = tracks->begin();
132 std::vector<Trk::VxTrackAtVertex *>::const_iterator ve = tracks->end();
133 unsigned int n_failed = 0;
134 std::vector<double> in_weights ( 0 );
135 std::vector<double> out_weights ( 0 );
136 std::vector<double> pu_weights ( 0 );
137 std::vector<double> no_correspondance ( 0 );
138
139
140 for ( ;
vt!=ve;++
vt ) {
141
142
143 if ( ( *vt ) !=nullptr ) {
144
145 ITrackLink * origLink = ( **vt ).trackOrParticleLink();
146
147 if ( origLink !=nullptr ) {
148
149 LinkToTrackParticleBase * tr_part = dynamic_cast< LinkToTrackParticleBase * > ( origLink );
150 if ( tr_part !=nullptr && tr_part->isValid()) {
151
152
153 std::map< Rec::TrackParticleTruthKey, TrackParticleTruth>::const_iterator ttItr = trackParticleTruthCollection->end();
154
155 const Rec::TrackParticle*
tp =
dynamic_cast<const Rec::TrackParticle*
>(**tr_part);
156 if(tp) {
157 const Rec::TrackParticleContainer* tpCont = dynamic_cast<const Rec::TrackParticleContainer*>(tr_part->getStorableObjectPointer());
158 if (tpCont) {
159 ElementLink<Rec::TrackParticleContainer> linkTruth;
162 ttItr = trackParticleTruthCollection->find(Rec::TrackParticleTruthKey(linkTruth));
163 }
164 }
165
166
167 if (ttItr != trackParticleTruthCollection->end() ) {
168 const HepMcParticleLink& particleLink = ttItr->second.particleLink();
169#ifdef HEPMC3
171#else
173#endif
174
175 if(genParticle) {
176 const auto *tpEvent = genParticle->parent_event();
177 if(tpEvent==genEvent) {
179 if (genParticle) pVertex = genParticle->production_vertex();
180 if ( pVertex) {
181 int link_pid = genParticle->pdg_id();
182 bool primary_track = false;
183 bool secondary_track = false;
184
185
186 do {
187#ifdef HEPMC3
188 auto idf_res = vertex_ids.find ( pVertex->id() );
189#else
191#endif
192
193
194
195 if ( idf_res != vertex_ids.end() ) {
196
197 primary_track = true;
198 in_weights.push_back ( ( *vt )->weight() );
199 }else {
200
201
202
203 if ( pVertex->particles_in_size() == 1 ) {
204
205#ifdef HEPMC3
206 auto inp = pVertex->particles_in()[0] ;
207#else
208 auto inp =*(pVertex->particles_in_const_begin()) ;
209#endif
210 auto tmpVertex_loc = inp ->production_vertex();
211 if ( inp ->pdg_id() == link_pid && tmpVertex_loc) {
212
213
214 pVertex = tmpVertex_loc;
215 }else {
216 secondary_track = true;
217 out_weights.push_back ( ( *vt )->weight() );
218 }
219 }else {
220
221
222
223 secondary_track = true;
224 out_weights.push_back ( ( *vt )->weight() );
225 }
226 }
227 }while ( !primary_track && !secondary_track );
228 }else {
229
230
231
233 msg() <<
"A truth particle with no production vertex found."<<
endmsg;
234 msg() <<
"Since it does not come from PV, consider as PileUp"<<
endmsg;
235 }
236 pu_weights.push_back ( ( *vt )->weight() );
237 }
238 }else {
239
240
241 pu_weights.push_back ( ( *vt )->weight() );
242 }
243 }else{
244
245
246
247 no_correspondance.push_back ( ( *vt )->weight() );
248 }
249 }
250 }else{
251 if (
msgLvl(MSG::DEBUG))
msg() <<
"This track at vertex has no valid link to its original trackparticle "<<
endmsg;
252 ++ n_failed;
253 }
254 }else{
255 if (
msgLvl(MSG::DEBUG))
msg() <<
"There is an empty pointer in the vector<VxTrackAtVertex *> for this vertex"<<
endmsg;
256 ++n_failed;
257 }
258 }
259 }
260
261
262
263
264
265
266
267
268
269 return new Trk::TrkPriVxPurity ( tracks->size(), pu_weights, no_correspondance,
270 n_failed, in_weights, out_weights );
271 }
272 }else{
274 msg()<<
"Empty vertex pointer received"<<
endmsg;
275 msg()<<
"Empty pointer returned."<<
endmsg;
276 }
277 return nullptr;
278 }
279 return nullptr;
280 }
Athena::TPCnvVers::Old TrackParticleTruthCollection
ServiceHandle< StoreGateSvc > & evtStore()
bool msgLvl(const MSG::Level lvl) const
DataModel_detail::const_iterator< DataVector > const_iterator
bool setElement(ElementType element)
Set to point to an element.
bool setStorableObject(BaseConstReference data, bool replace=false, IProxyDict *sg=0)
Set link to point to a new container (storable).
HepMC::ConstGenParticlePtr cptr() const
Dereference.
HepMC::ConstGenParticlePtr scptr() const
Dereference/smart pointer.
::StatusCode StatusCode
StatusCode definition for legacy code.
const GenParticle * ConstGenParticlePtr
const HepMC::GenVertex * ConstGenVertexPtr