63 ATH_MSG_FATAL(
"Incorrect number of Psi daughters (should be 3 or 4)");
64 return StatusCode::FAILURE;
67 constexpr int topoN = 3;
70 return StatusCode::FAILURE;
72 std::array<SG::WriteHandle<xAOD::VertexContainer>, topoN> VtxWriteHandles;
int ikey(0);
75 ATH_CHECK( VtxWriteHandles[ikey].record(std::make_unique<xAOD::VertexContainer>(), std::make_unique<xAOD::VertexAuxContainer>()) );
84 if (pvContainer.
cptr()->size()==0) {
86 return StatusCode::RECOVERABLE;
95 ATH_CHECK( refPvContainer.
record(std::make_unique<xAOD::VertexContainer>(), std::make_unique<xAOD::VertexAuxContainer>()) );
98 std::vector<Trk::VxCascadeInfo*> cascadeinfoContainer;
99 std::vector<Trk::VxCascadeInfo*> cascadeinfoContainer_noConstr;
108 SG::AuxElement::Decorator<VertexLinkVector> CascadeLinksDecor(
"CascadeVertexLinks");
109 SG::AuxElement::Decorator<VertexLinkVector> Psi1LinksDecor(
"Psi1VertexLinks");
110 SG::AuxElement::Decorator<VertexLinkVector> Psi2LinksDecor(
"Psi2VertexLinks");
111 SG::AuxElement::Decorator<float> chi2_decor(
"ChiSquared");
112 SG::AuxElement::Decorator<int> ndof_decor(
"nDoF");
113 SG::AuxElement::Decorator<float> chi2_nc_decor(
"ChiSquared_nc");
114 SG::AuxElement::Decorator<int> ndof_nc_decor(
"nDoF_nc");
115 SG::AuxElement::Decorator<float> Pt_decor(
"Pt");
116 SG::AuxElement::Decorator<float> PtErr_decor(
"PtErr");
117 SG::AuxElement::Decorator<float> chi2_SV1_decor(
"ChiSquared_SV1");
118 SG::AuxElement::Decorator<float> chi2_nc_SV1_decor(
"ChiSquared_nc_SV1");
119 SG::AuxElement::Decorator<float> chi2_V1_decor(
"ChiSquared_V1");
120 SG::AuxElement::Decorator<int> ndof_V1_decor(
"nDoF_V1");
121 SG::AuxElement::Decorator<float> lxy_SV1_decor(
"lxy_SV1");
122 SG::AuxElement::Decorator<float> lxyErr_SV1_decor(
"lxyErr_SV1");
123 SG::AuxElement::Decorator<float> a0xy_SV1_decor(
"a0xy_SV1");
124 SG::AuxElement::Decorator<float> a0xyErr_SV1_decor(
"a0xyErr_SV1");
125 SG::AuxElement::Decorator<float> a0z_SV1_decor(
"a0z_SV1");
126 SG::AuxElement::Decorator<float> a0zErr_SV1_decor(
"a0zErr_SV1");
127 SG::AuxElement::Decorator<float> chi2_SV2_decor(
"ChiSquared_SV2");
128 SG::AuxElement::Decorator<float> chi2_nc_SV2_decor(
"ChiSquared_nc_SV2");
129 SG::AuxElement::Decorator<float> chi2_V2_decor(
"ChiSquared_V2");
130 SG::AuxElement::Decorator<int> ndof_V2_decor(
"nDoF_V2");
131 SG::AuxElement::Decorator<float> lxy_SV2_decor(
"lxy_SV2");
132 SG::AuxElement::Decorator<float> lxyErr_SV2_decor(
"lxyErr_SV2");
133 SG::AuxElement::Decorator<float> a0xy_SV2_decor(
"a0xy_SV2");
134 SG::AuxElement::Decorator<float> a0xyErr_SV2_decor(
"a0xyErr_SV2");
135 SG::AuxElement::Decorator<float> a0z_SV2_decor(
"a0z_SV2");
136 SG::AuxElement::Decorator<float> a0zErr_SV2_decor(
"a0zErr_SV2");
144 for(
size_t ic=0; ic<cascadeinfoContainer.size(); ic++) {
146 if(cascade_info==
nullptr) {
153 const std::vector<xAOD::Vertex*> &cascadeVertices = cascade_info->
vertices();
154 if(cascadeVertices.size() != topoN)
ATH_MSG_ERROR(
"Incorrect number of vertices");
155 if(cascadeVertices[0]==
nullptr || cascadeVertices[1]==
nullptr || cascadeVertices[2]==
nullptr)
ATH_MSG_ERROR(
"Error null vertex");
157 for(
int i=0; i<topoN; i++) VtxWriteHandles[i].ptr()->push_back(cascadeVertices[i]);
160 const auto mainVertex = cascadeVertices[2];
161 const std::vector< std::vector<TLorentzVector> > &moms = cascade_info->
getParticleMoms();
164 std::vector<VertexLink> precedingVertexLinks;
168 if( vertexLink1.
isValid() ) precedingVertexLinks.push_back( vertexLink1 );
172 if( vertexLink2.
isValid() ) precedingVertexLinks.push_back( vertexLink2 );
173 CascadeLinksDecor(*mainVertex) = precedingVertexLinks;
185 std::vector<const xAOD::Vertex*> psi2VerticestoLink;
186 if(psi2Vertex) psi2VerticestoLink.push_back(psi2Vertex);
190 std::vector<const xAOD::Vertex*> psi1VerticestoLink;
191 if(psi1Vertex) psi1VerticestoLink.push_back(psi1Vertex);
212 chi2_decor(*mainVertex) = cascade_info->
fitChi2();
213 ndof_decor(*mainVertex) = cascade_info->
nDoF();
214 chi2_nc_decor(*mainVertex) = cascade_info_noConstr ? cascade_info_noConstr->
fitChi2() : -999999.;
215 ndof_nc_decor(*mainVertex) = cascade_info_noConstr ? cascade_info_noConstr->
nDoF() : -1;
218 chi2_SV1_decor(*cascadeVertices[0]) =
m_V0Tools->chisq(cascadeVertices[0]);
219 chi2_nc_SV1_decor(*cascadeVertices[0]) = cascade_info_noConstr ?
m_V0Tools->chisq(cascade_info_noConstr->
vertices()[0]) : -999999.;
220 chi2_V1_decor(*cascadeVertices[0]) =
m_V0Tools->chisq(psi1Vertex);
221 ndof_V1_decor(*cascadeVertices[0]) =
m_V0Tools->ndof(psi1Vertex);
222 lxy_SV1_decor(*cascadeVertices[0]) =
m_CascadeTools->lxy(moms[0],cascadeVertices[0],mainVertex);
223 lxyErr_SV1_decor(*cascadeVertices[0]) =
m_CascadeTools->lxyError(moms[0],cascade_info->
getCovariance()[0],cascadeVertices[0],mainVertex);
224 a0z_SV1_decor(*cascadeVertices[0]) =
m_CascadeTools->a0z(moms[0],cascadeVertices[0],mainVertex);
225 a0zErr_SV1_decor(*cascadeVertices[0]) =
m_CascadeTools->a0zError(moms[0],cascade_info->
getCovariance()[0],cascadeVertices[0],mainVertex);
226 a0xy_SV1_decor(*cascadeVertices[0]) =
m_CascadeTools->a0xy(moms[0],cascadeVertices[0],mainVertex);
227 a0xyErr_SV1_decor(*cascadeVertices[0]) =
m_CascadeTools->a0xyError(moms[0],cascade_info->
getCovariance()[0],cascadeVertices[0],mainVertex);
230 chi2_SV2_decor(*cascadeVertices[1]) =
m_V0Tools->chisq(cascadeVertices[1]);
231 chi2_nc_SV2_decor(*cascadeVertices[1]) = cascade_info_noConstr ?
m_V0Tools->chisq(cascade_info_noConstr->
vertices()[1]) : -999999.;
232 chi2_V2_decor(*cascadeVertices[1]) =
m_V0Tools->chisq(psi2Vertex);
233 ndof_V2_decor(*cascadeVertices[1]) =
m_V0Tools->ndof(psi2Vertex);
234 lxy_SV2_decor(*cascadeVertices[1]) =
m_CascadeTools->lxy(moms[1],cascadeVertices[1],mainVertex);
235 lxyErr_SV2_decor(*cascadeVertices[1]) =
m_CascadeTools->lxyError(moms[1],cascade_info->
getCovariance()[1],cascadeVertices[1],mainVertex);
236 a0z_SV2_decor(*cascadeVertices[1]) =
m_CascadeTools->a0z(moms[1],cascadeVertices[1],mainVertex);
237 a0zErr_SV2_decor(*cascadeVertices[1]) =
m_CascadeTools->a0zError(moms[1],cascade_info->
getCovariance()[1],cascadeVertices[1],mainVertex);
238 a0xy_SV2_decor(*cascadeVertices[1]) =
m_CascadeTools->a0xy(moms[1],cascadeVertices[1],mainVertex);
239 a0xyErr_SV2_decor(*cascadeVertices[1]) =
m_CascadeTools->a0xyError(moms[1],cascade_info->
getCovariance()[1],cascadeVertices[1],mainVertex);
247 for (
auto cascade_info : cascadeinfoContainer)
delete cascade_info;
248 for (
auto cascade_info_noConstr : cascadeinfoContainer_noConstr)
delete cascade_info_noConstr;
250 return StatusCode::SUCCESS;
366 StatusCode
PsiPlusPsiCascade::performSearch(std::vector<Trk::VxCascadeInfo*> *cascadeinfoContainer, std::vector<Trk::VxCascadeInfo*> *cascadeinfoContainer_noConstr,
const EventContext& ctx)
const {
368 assert(cascadeinfoContainer!=
nullptr && cascadeinfoContainer_noConstr!=
nullptr);
374 std::vector<const xAOD::TrackParticle*> tracksJpsi1;
375 std::vector<const xAOD::TrackParticle*> tracksJpsi2;
376 std::vector<const xAOD::TrackParticle*> tracksDiTrk1;
377 std::vector<const xAOD::TrackParticle*> tracksDiTrk2;
378 std::vector<const xAOD::TrackParticle*> tracksPsi1;
379 std::vector<const xAOD::TrackParticle*> tracksPsi2;
380 std::vector<double> massesPsi1;
385 std::vector<double> massesPsi2;
400 std::vector<const xAOD::Vertex*> selectedPsi2Candidates;
401 for(
auto vxcItr=psi2Container.
cptr()->cbegin(); vxcItr!=psi2Container.
cptr()->cend(); ++vxcItr) {
407 if(flagAcc.isAvailable(*vtx) && flagAcc(*vtx)) {
414 double mass_psi2 =
m_V0Tools->invariantMass(*vxcItr, massesPsi2);
415 if (mass_psi2 < m_psi2MassLower || mass_psi2 >
m_psi2MassUpper)
continue;
418 TLorentzVector p4_mu1, p4_mu2;
419 p4_mu1.SetPtEtaPhiM( (*vxcItr)->trackParticle(0)->pt(),
420 (*vxcItr)->trackParticle(0)->eta(),
422 p4_mu2.SetPtEtaPhiM( (*vxcItr)->trackParticle(1)->pt(),
423 (*vxcItr)->trackParticle(1)->eta(),
425 double mass_jpsi2 = (p4_mu1 + p4_mu2).M();
426 if (mass_jpsi2 < m_jpsi2MassLower || mass_jpsi2 >
m_jpsi2MassUpper)
continue;
429 TLorentzVector p4_trk1, p4_trk2;
430 p4_trk1.SetPtEtaPhiM( (*vxcItr)->trackParticle(2)->pt(),
431 (*vxcItr)->trackParticle(2)->eta(),
433 p4_trk2.SetPtEtaPhiM( (*vxcItr)->trackParticle(3)->pt(),
434 (*vxcItr)->trackParticle(3)->eta(),
436 double mass_diTrk2 = (p4_trk1 + p4_trk2).M();
440 double chi2DOF = (*vxcItr)->chiSquared()/(*vxcItr)->numberDoF();
443 selectedPsi2Candidates.push_back(*vxcItr);
445 if(selectedPsi2Candidates.size()==0)
return StatusCode::SUCCESS;
448 std::vector<const xAOD::Vertex*> selectedPsi1Candidates;
449 for(
auto vxcItr=psi1Container.
cptr()->cbegin(); vxcItr!=psi1Container.
cptr()->cend(); ++vxcItr) {
455 if(flagAcc.isAvailable(*vtx) && flagAcc(*vtx)) {
462 double mass_psi1 =
m_V0Tools->invariantMass(*vxcItr,massesPsi1);
463 if(mass_psi1 < m_psi1MassLower || mass_psi1 >
m_psi1MassUpper)
continue;
466 TLorentzVector p4_mu1, p4_mu2;
467 p4_mu1.SetPtEtaPhiM( (*vxcItr)->trackParticle(0)->pt(),
468 (*vxcItr)->trackParticle(0)->eta(),
470 p4_mu2.SetPtEtaPhiM( (*vxcItr)->trackParticle(1)->pt(),
471 (*vxcItr)->trackParticle(1)->eta(),
473 double mass_jpsi1 = (p4_mu1 + p4_mu2).M();
474 if (mass_jpsi1 < m_jpsi1MassLower || mass_jpsi1 >
m_jpsi1MassUpper)
continue;
477 TLorentzVector p4_trk1, p4_trk2;
478 p4_trk1.SetPtEtaPhiM( (*vxcItr)->trackParticle(2)->pt(),
479 (*vxcItr)->trackParticle(2)->eta(),
481 p4_trk2.SetPtEtaPhiM( (*vxcItr)->trackParticle(3)->pt(),
482 (*vxcItr)->trackParticle(3)->eta(),
484 double mass_diTrk1 = (p4_trk1 + p4_trk2).M();
488 double chi2DOF = (*vxcItr)->chiSquared()/(*vxcItr)->numberDoF();
491 selectedPsi1Candidates.push_back(*vxcItr);
493 if(selectedPsi1Candidates.size()==0)
return StatusCode::SUCCESS;
495 std::vector<std::pair<const xAOD::Vertex*, const xAOD::Vertex*> > candidatePairs;
496 for(
auto psi1Itr=selectedPsi1Candidates.cbegin(); psi1Itr!=selectedPsi1Candidates.cend(); ++psi1Itr) {
498 for(
size_t i=0; i<(*psi1Itr)->nTrackParticles(); i++) tracksPsi1.push_back((*psi1Itr)->trackParticle(i));
499 for(
auto psi2Itr=selectedPsi2Candidates.cbegin(); psi2Itr!=selectedPsi2Candidates.cend(); ++psi2Itr) {
501 for(
size_t j=0; j<(*psi2Itr)->nTrackParticles(); j++) {
502 if(std::find(tracksPsi1.cbegin(), tracksPsi1.cend(), (*psi2Itr)->trackParticle(j)) != tracksPsi1.cend()) {
skip =
true;
break; }
506 for(
size_t ic=0; ic<candidatePairs.size(); ic++) {
507 const xAOD::Vertex* psi1Vertex = candidatePairs[ic].first;
508 const xAOD::Vertex* psi2Vertex = candidatePairs[ic].second;
509 if((psi1Vertex == *psi1Itr && psi2Vertex == *psi2Itr) || (psi1Vertex == *psi2Itr && psi2Vertex == *psi1Itr)) {
skip =
true;
break; }
513 candidatePairs.push_back(std::pair<const xAOD::Vertex*, const xAOD::Vertex*>(*psi1Itr,*psi2Itr));
517 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(); } );
519 candidatePairs.erase(candidatePairs.begin()+
m_maxCandidates, candidatePairs.end());
522 for(
size_t ic=0; ic<candidatePairs.size(); ic++) {
523 const xAOD::Vertex* psi1Vertex = candidatePairs[ic].first;
524 const xAOD::Vertex* psi2Vertex = candidatePairs[ic].second;
528 if (tracksPsi1.size() != massesPsi1.size()) {
529 ATH_MSG_ERROR(
"Problems with Psi1 input: number of tracks or track mass inputs is not correct!");
533 if (tracksPsi2.size() != massesPsi2.size()) {
534 ATH_MSG_ERROR(
"Problems with Psi2 input: number of tracks or track mass inputs is not correct!");
540 tracksDiTrk1.clear();
548 tracksDiTrk2.clear();
554 TLorentzVector p4_moth;
567 std::unique_ptr<Trk::IVKalState> state =
m_iVertexFitter->makeState(ctx);
573 std::vector<Trk::VertexID> vrtList;
582 vrtList.push_back(vID1);
590 vrtList.push_back(vID2);
592 std::vector<const xAOD::TrackParticle*> tp; tp.clear();
593 std::vector<double> tp_masses; tp_masses.clear();
596 std::vector<Trk::VertexID> cnstV; cnstV.clear();
602 std::vector<Trk::VertexID> cnstV; cnstV.clear();
608 std::vector<Trk::VertexID> cnstV; cnstV.clear();
614 std::vector<Trk::VertexID> cnstV; cnstV.clear();
620 std::unique_ptr<Trk::VxCascadeInfo> result(
m_iVertexFitter->fitCascade(*state));
623 if (result !=
nullptr) {
624 for(
auto v : result->vertices()) {
625 if(v->nTrackParticles()==0) {
626 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
627 v->setTrackParticleLinks(nullLinkVector);
634 result->setSVOwnership(
true);
637 double chi2DOF = result->fitChi2()/result->nDoF();
641 cascadeinfoContainer->push_back(result.release());
649 std::unique_ptr<Trk::IVKalState> state (
m_iVertexFitter->makeState(ctx));
651 std::vector<Trk::VertexID> vrtList_nc;
654 vrtList_nc.push_back(vID1_nc);
657 vrtList_nc.push_back(vID2_nc);
659 std::vector<const xAOD::TrackParticle*> tp; tp.clear();
660 std::vector<double> tp_masses; tp_masses.clear();
663 std::unique_ptr<Trk::VxCascadeInfo> result_nc(
m_iVertexFitter->fitCascade(*state));
665 if (result_nc !=
nullptr) {
666 for(
auto v : result_nc->vertices()) {
667 if(v->nTrackParticles()==0) {
668 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
669 v->setTrackParticleLinks(nullLinkVector);
676 result_nc->setSVOwnership(
true);
677 cascadeinfoContainer_noConstr->push_back(result_nc.release());
679 else cascadeinfoContainer_noConstr->push_back(0);
681 else cascadeinfoContainer_noConstr->push_back(0);
685 return StatusCode::SUCCESS;