167 {
170 return StatusCode::FAILURE;
171 }
172
173 constexpr int topoN = 3;
176 return StatusCode::FAILURE;
177 }
178 std::array<SG::WriteHandle<xAOD::VertexContainer>, topoN> VtxWriteHandles; int ikey(0);
179 for(
const SG::WriteHandleKey<xAOD::VertexContainer>& key :
m_outputsKeys) {
180 VtxWriteHandles[ikey] = SG::WriteHandle<xAOD::VertexContainer>(key, ctx);
181 ATH_CHECK( VtxWriteHandles[ikey].record(std::make_unique<xAOD::VertexContainer>(), std::make_unique<xAOD::VertexAuxContainer>()) );
182 ikey++;
183 }
184
185
186
187
190 if (pvContainer.cptr()->size()==0) {
191 ATH_MSG_WARNING(
"You have no primary vertices: " << pvContainer.cptr()->size());
192 return StatusCode::RECOVERABLE;
193 }
194
195
198
199 std::vector<const xAOD::TrackParticle*> tracksJpsi1;
200 std::vector<const xAOD::TrackParticle*> tracksDiTrk1;
201 std::vector<const xAOD::TrackParticle*> tracksPsi1;
202 std::vector<const xAOD::TrackParticle*> tracksJpsi2;
203 std::vector<const xAOD::TrackParticle*> tracksDiTrk2;
204 std::vector<const xAOD::TrackParticle*> tracksPsi2;
205 std::vector<const xAOD::TrackParticle*> inputTracks;
206 std::vector<double> massesPsi1;
211 std::vector<double> massesPsi2;
216 std::vector<double> massesInputTracks;
217 massesInputTracks.reserve(massesPsi1.size() + massesPsi2.size());
218 for(auto mass : massesPsi1) massesInputTracks.push_back(mass);
219 for(auto mass : massesPsi2) massesInputTracks.push_back(mass);
220
221
226
227
228 std::vector<const xAOD::Vertex*> selectedPsi1Candidates;
229 for(auto vxcItr=psi1Container->cbegin(); vxcItr!=psi1Container->cend(); ++vxcItr) {
230
235 if(flagAcc.isAvailable(*vtx) && flagAcc(*vtx)) {
237 }
238 }
240
241
243 double mass_psi1 =
m_V0Tools->invariantMass(*vxcItr, massesPsi1);
244 if (mass_psi1 < m_psi1MassLower || mass_psi1 >
m_psi1MassUpper)
continue;
245 }
246
247 TLorentzVector p4_mu1, p4_mu2;
248 p4_mu1.SetPtEtaPhiM( (*vxcItr)->trackParticle(0)->pt(),
249 (*vxcItr)->trackParticle(0)->eta(),
251 p4_mu2.SetPtEtaPhiM( (*vxcItr)->trackParticle(1)->pt(),
252 (*vxcItr)->trackParticle(1)->eta(),
254 double mass_jpsi1 = (p4_mu1 + p4_mu2).M();
255 if (mass_jpsi1 < m_jpsi1MassLower || mass_jpsi1 >
m_jpsi1MassUpper)
continue;
256
258 TLorentzVector p4_trk1, p4_trk2;
259 p4_trk1.SetPtEtaPhiM( (*vxcItr)->trackParticle(2)->pt(),
260 (*vxcItr)->trackParticle(2)->eta(),
262 p4_trk2.SetPtEtaPhiM( (*vxcItr)->trackParticle(3)->pt(),
263 (*vxcItr)->trackParticle(3)->eta(),
265 double mass_diTrk1 = (p4_trk1 + p4_trk2).M();
267 }
268
269 selectedPsi1Candidates.push_back(*vxcItr);
270 }
271
272
273 std::vector<const xAOD::Vertex*> selectedPsi2Candidates;
274 for(auto vxcItr=psi2Container->cbegin(); vxcItr!=psi2Container->cend(); ++vxcItr) {
275
280 if(flagAcc.isAvailable(*vtx) && flagAcc(*vtx)) {
282 }
283 }
285
286
288 double mass_psi2 =
m_V0Tools->invariantMass(*vxcItr,massesPsi2);
289 if(mass_psi2 < m_psi2MassLower || mass_psi2 >
m_psi2MassUpper)
continue;
290 }
291
292 TLorentzVector p4_mu1, p4_mu2;
293 p4_mu1.SetPtEtaPhiM( (*vxcItr)->trackParticle(0)->pt(),
294 (*vxcItr)->trackParticle(0)->eta(),
296 p4_mu2.SetPtEtaPhiM( (*vxcItr)->trackParticle(1)->pt(),
297 (*vxcItr)->trackParticle(1)->eta(),
299 double mass_jpsi2 = (p4_mu1 + p4_mu2).M();
300 if (mass_jpsi2 < m_jpsi2MassLower || mass_jpsi2 >
m_jpsi2MassUpper)
continue;
301
303 TLorentzVector p4_trk1, p4_trk2;
304 p4_trk1.SetPtEtaPhiM( (*vxcItr)->trackParticle(2)->pt(),
305 (*vxcItr)->trackParticle(2)->eta(),
307 p4_trk2.SetPtEtaPhiM( (*vxcItr)->trackParticle(3)->pt(),
308 (*vxcItr)->trackParticle(3)->eta(),
310 double mass_diTrk2 = (p4_trk1 + p4_trk2).M();
312 }
313 selectedPsi2Candidates.push_back(*vxcItr);
314 }
315
316 std::vector<std::pair<const xAOD::Vertex*, const xAOD::Vertex*> > candidatePairs;
317 for(auto psi1Itr=selectedPsi1Candidates.cbegin(); psi1Itr!=selectedPsi1Candidates.cend(); ++psi1Itr) {
318 tracksPsi1.clear();
319 tracksPsi1.reserve((*psi1Itr)->nTrackParticles());
320 for(
size_t i=0;
i<(*psi1Itr)->nTrackParticles();
i++) tracksPsi1.push_back((*psi1Itr)->trackParticle(i));
321 for(auto psi2Itr=selectedPsi2Candidates.cbegin(); psi2Itr!=selectedPsi2Candidates.cend(); ++psi2Itr) {
323 for(
size_t j=0;
j<(*psi2Itr)->nTrackParticles();
j++) {
324 if(std::find(tracksPsi1.cbegin(), tracksPsi1.cend(), (*psi2Itr)->trackParticle(j)) != tracksPsi1.cend()) {
skip =
true;
break; }
325 }
328 for(
size_t ic=0;
ic<candidatePairs.size();
ic++) {
331 if((psi1Vertex == *psi1Itr && psi2Vertex == *psi2Itr) || (psi1Vertex == *psi2Itr && psi2Vertex == *psi1Itr)) {
skip =
true;
break; }
332 }
333 }
335 candidatePairs.push_back(std::pair<const xAOD::Vertex*, const xAOD::Vertex*>(*psi1Itr,*psi2Itr));
336 }
337 }
338
339 std::sort( candidatePairs.begin(), candidatePairs.end(), [](std::pair<const xAOD::Vertex*, const xAOD::Vertex*>
a, std::pair<const xAOD::Vertex*, const xAOD::Vertex*> b) { return a.first->chiSquared()/a.first->numberDoF()+a.second->chiSquared()/a.second->numberDoF() < b.first->chiSquared()/b.first->numberDoF()+b.second->chiSquared()/b.second->numberDoF(); } );
341 candidatePairs.erase(candidatePairs.begin()+
m_maxCandidates, candidatePairs.end());
342 }
343
344 for(
size_t ic=0;
ic<candidatePairs.size();
ic++) {
347
348 tracksPsi1.clear();
351 if (tracksPsi1.size() != massesPsi1.size()) {
352 ATH_MSG_ERROR(
"Problems with Psi1 input: number of tracks or track mass inputs is not correct!");
353 }
354 tracksPsi2.clear();
357 if (tracksPsi2.size() != massesPsi2.size()) {
358 ATH_MSG_ERROR(
"Problems with Psi2 input: number of tracks or track mass inputs is not correct!");
359 }
360
361 tracksJpsi1.clear();
364 tracksDiTrk1.clear();
368 }
369 tracksJpsi2.clear();
372 tracksDiTrk2.clear();
376 }
377
378 TLorentzVector p4_moth;
383 }
387 }
389
390 inputTracks.clear();
399
400
401 std::unique_ptr<Trk::IVKalState> state =
m_iVertexFitter->makeState(ctx);
405 }
408 }
411 }
414 }
417 }
420 }
421
422
423 Amg::Vector3D startingPoint((psi1Vertex->
x()+psi2Vertex->
x())/2,(psi1Vertex->
y()+psi2Vertex->
y())/2,(psi1Vertex->
z()+psi2Vertex->
z())/2);
424
425
426 std::unique_ptr<xAOD::Vertex> theResult(
m_iVertexFitter->fit(inputTracks, startingPoint, *state) );
427
428 if(theResult != nullptr){
429
430 double chi2DOF = theResult->chiSquared()/theResult->numberDoF();
432 if(chi2CutPassed) {
434 for(
size_t i=0;
i<theResult->trackParticleLinks().
size();
i++) {
436 mylink.setStorableObject(*trackContainer.ptr(), true);
437 tpLinkVector.push_back( mylink );
438 }
439 theResult->clearTracks();
440 theResult->setTrackParticleLinks( tpLinkVector );
441
443 xAOD::BPhysHelper psi1_helper(psi1Vertex_);
444 psi1_helper.setRefTrks();
448 mylink.setStorableObject(*trackContainer.ptr(), true);
449 tpLinkVector_psi1.push_back( mylink );
450 }
453
455 xAOD::BPhysHelper psi2_helper(psi2Vertex_);
456 psi2_helper.setRefTrks();
460 mylink.setStorableObject(*trackContainer.ptr(), true);
461 tpLinkVector_psi2.push_back( mylink );
462 }
465
466 VtxWriteHandles[0].ptr()->push_back(psi1Vertex_);
467 VtxWriteHandles[1].ptr()->push_back(psi2Vertex_);
468
469
474 if( vertexLink1.
isValid() ) precedingVertexLinks.push_back( vertexLink1 );
478 if( vertexLink2.
isValid() ) precedingVertexLinks.push_back( vertexLink2 );
479
480 SG::AuxElement::Decorator<VertexLinkVector> PrecedingLinksDecor("PrecedingVertexLinks");
481 PrecedingLinksDecor(*theResult.get()) = precedingVertexLinks;
482
483 xAOD::BPhysHypoHelper vtx(
m_hypoName, theResult.get());
484 vtx.setRefTrks();
485 vtx.setPass(true);
486
487 VtxWriteHandles[2].ptr()->push_back(theResult.release());
488 }
489 }
490 }
491
496
499 ATH_CHECK( refPvContainer.record(std::make_unique<xAOD::VertexContainer>(), std::make_unique<xAOD::VertexAuxContainer>()) );
500
501 if(VtxWriteHandles[2]->
size()>0) {
502 for(
int i=0;
i<topoN;
i++) {
505 ATH_MSG_FATAL(
"FillCandwithRefittedVertices failed - check the vertices you passed");
507 }
508 }
509 }
510 }
511 else {
512 if(VtxWriteHandles[2]->
size()>0) {
513 for(
int i=0;
i<topoN;
i++) {
516 ATH_MSG_FATAL(
"FillCandExistingVertices failed - check the vertices you passed");
518 }
519 }
520 }
521 }
522
523 return StatusCode::SUCCESS;
524 }
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_WARNING(x)
size_t size() const
Number of registered mappings.
bool setElement(ElementType element)
Set link to point to an Element (slowest).
bool isValid() const
Check if the element can be found.
bool setStorableObject(BaseConstReference data, bool replace=false)
Set link storable to data object pointed by data (slower).
virtual double phi() const override final
The azimuthal angle ( ) of the particle (has range to .).
virtual double pt() const override final
The transverse momentum ( ) of the particle.
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
float z() const
Returns the z position.
void setTrackParticleLinks(const TrackParticleLinks_t &trackParticles)
Set all track particle links at once.
size_t nTrackParticles() const
Get the number of tracks associated with this vertex.
void clearTracks()
Remove all tracks from the vertex.
const TrackParticleLinks_t & trackParticleLinks() const
Get all the particles associated with the vertex.
const TrackParticle * trackParticle(size_t i) const
Get the pointer to a given track that was used in vertex reco.
float y() const
Returns the y position.
float x() const
Returns the x position.
Eigen::Matrix< double, 3, 1 > Vector3D
ElementLink< xAOD::TrackParticleContainer > TrackParticleLink
ElementLink< xAOD::VertexContainer > VertexLink
std::vector< VertexLink > VertexLinkVector
std::vector< TrackParticleLink > TrackParticleLinkVector
::StatusCode StatusCode
StatusCode definition for legacy code.
const float SC[NC]
Cross sections for Carbon.
float j(const xAOD::IParticle &, const xAOD::TrackMeasurementValidation &hit, const Eigen::Matrix3d &jab_inv)
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
Vertex_v1 Vertex
Define the latest version of the vertex class.