10 #include "GaudiKernel/IPartPropSvc.h"
15 #include "HepPDT/ParticleDataTable.hh"
45 IPartPropSvc* partPropSvc =
nullptr;
46 ATH_CHECK( service(
"PartPropSvc", partPropSvc,
true) );
47 auto pdt = partPropSvc->PDT();
65 return StatusCode::SUCCESS;
70 ATH_MSG_FATAL(
"Incorrect number of Psi daughters (should be 3 or 4)");
71 return StatusCode::FAILURE;
74 constexpr
int topoN = 3;
77 return StatusCode::FAILURE;
79 std::array<SG::WriteHandle<xAOD::VertexContainer>, topoN> VtxWriteHandles;
int ikey(0);
82 ATH_CHECK( VtxWriteHandles[ikey].record(std::make_unique<xAOD::VertexContainer>(), std::make_unique<xAOD::VertexAuxContainer>()) );
93 return StatusCode::RECOVERABLE;
102 ATH_CHECK( refPvContainer.
record(std::make_unique<xAOD::VertexContainer>(), std::make_unique<xAOD::VertexAuxContainer>()) );
105 std::vector<Trk::VxCascadeInfo*> cascadeinfoContainer;
106 std::vector<Trk::VxCascadeInfo*> cascadeinfoContainer_noConstr;
151 for(
size_t ic=0;
ic<cascadeinfoContainer.size();
ic++) {
153 if(cascade_info==
nullptr)
ATH_MSG_ERROR(
"CascadeInfo is null");
157 const std::vector<xAOD::Vertex*> &cascadeVertices = cascade_info->
vertices();
158 if(cascadeVertices.size() != topoN)
ATH_MSG_ERROR(
"Incorrect number of vertices");
159 if(cascadeVertices[0]==
nullptr || cascadeVertices[1]==
nullptr || cascadeVertices[2]==
nullptr)
ATH_MSG_ERROR(
"Error null vertex");
161 for(
int i=0;
i<topoN;
i++) VtxWriteHandles[
i].ptr()->push_back(cascadeVertices[
i]);
164 const auto mainVertex = cascadeVertices[2];
165 const std::vector< std::vector<TLorentzVector> > &moms = cascade_info->
getParticleMoms();
168 std::vector<VertexLink> precedingVertexLinks;
172 if( vertexLink1.
isValid() ) precedingVertexLinks.push_back( vertexLink1 );
176 if( vertexLink2.
isValid() ) precedingVertexLinks.push_back( vertexLink2 );
177 CascadeLinksDecor(*mainVertex) = precedingVertexLinks;
181 if(
m_vtx2Daug_num==4) psi2Vertex = BPhysPVCascadeTools::FindVertex<4>(psi2Container.
cptr(), cascadeVertices[1]);
182 else psi2Vertex = BPhysPVCascadeTools::FindVertex<3>(psi2Container.
cptr(), cascadeVertices[1]);
185 if(
m_vtx1Daug_num==4) psi1Vertex = BPhysPVCascadeTools::FindVertex<4>(psi1Container.
cptr(), cascadeVertices[0]);
186 else psi1Vertex = BPhysPVCascadeTools::FindVertex<3>(psi1Container.
cptr(), cascadeVertices[0]);
189 std::vector<const xAOD::Vertex*> psi2VerticestoLink;
190 if(psi2Vertex) psi2VerticestoLink.push_back(psi2Vertex);
194 std::vector<const xAOD::Vertex*> psi1VerticestoLink;
195 if(psi1Vertex) psi1VerticestoLink.push_back(psi1Vertex);
216 chi2_decor(*mainVertex) = cascade_info->
fitChi2();
217 ndof_decor(*mainVertex) = cascade_info->
nDoF();
218 chi2_nc_decor(*mainVertex) = cascade_info_noConstr ? cascade_info_noConstr->
fitChi2() : -999999.;
219 ndof_nc_decor(*mainVertex) = cascade_info_noConstr ? cascade_info_noConstr->
nDoF() : -1;
222 chi2_SV1_decor(*cascadeVertices[0]) =
m_V0Tools->chisq(cascadeVertices[0]);
223 chi2_nc_SV1_decor(*cascadeVertices[0]) = cascade_info_noConstr ?
m_V0Tools->chisq(cascade_info_noConstr->
vertices()[0]) : -999999.;
224 chi2_V1_decor(*cascadeVertices[0]) =
m_V0Tools->chisq(psi1Vertex);
225 ndof_V1_decor(*cascadeVertices[0]) =
m_V0Tools->ndof(psi1Vertex);
226 lxy_SV1_decor(*cascadeVertices[0]) =
m_CascadeTools->lxy(moms[0],cascadeVertices[0],mainVertex);
227 lxyErr_SV1_decor(*cascadeVertices[0]) =
m_CascadeTools->lxyError(moms[0],cascade_info->
getCovariance()[0],cascadeVertices[0],mainVertex);
228 a0z_SV1_decor(*cascadeVertices[0]) =
m_CascadeTools->a0z(moms[0],cascadeVertices[0],mainVertex);
229 a0zErr_SV1_decor(*cascadeVertices[0]) =
m_CascadeTools->a0zError(moms[0],cascade_info->
getCovariance()[0],cascadeVertices[0],mainVertex);
230 a0xy_SV1_decor(*cascadeVertices[0]) =
m_CascadeTools->a0xy(moms[0],cascadeVertices[0],mainVertex);
231 a0xyErr_SV1_decor(*cascadeVertices[0]) =
m_CascadeTools->a0xyError(moms[0],cascade_info->
getCovariance()[0],cascadeVertices[0],mainVertex);
234 chi2_SV2_decor(*cascadeVertices[1]) =
m_V0Tools->chisq(cascadeVertices[1]);
235 chi2_nc_SV2_decor(*cascadeVertices[1]) = cascade_info_noConstr ?
m_V0Tools->chisq(cascade_info_noConstr->
vertices()[1]) : -999999.;
236 chi2_V2_decor(*cascadeVertices[1]) =
m_V0Tools->chisq(psi2Vertex);
237 ndof_V2_decor(*cascadeVertices[1]) =
m_V0Tools->ndof(psi2Vertex);
238 lxy_SV2_decor(*cascadeVertices[1]) =
m_CascadeTools->lxy(moms[1],cascadeVertices[1],mainVertex);
239 lxyErr_SV2_decor(*cascadeVertices[1]) =
m_CascadeTools->lxyError(moms[1],cascade_info->
getCovariance()[1],cascadeVertices[1],mainVertex);
240 a0z_SV2_decor(*cascadeVertices[1]) =
m_CascadeTools->a0z(moms[1],cascadeVertices[1],mainVertex);
241 a0zErr_SV2_decor(*cascadeVertices[1]) =
m_CascadeTools->a0zError(moms[1],cascade_info->
getCovariance()[1],cascadeVertices[1],mainVertex);
242 a0xy_SV2_decor(*cascadeVertices[1]) =
m_CascadeTools->a0xy(moms[1],cascadeVertices[1],mainVertex);
243 a0xyErr_SV2_decor(*cascadeVertices[1]) =
m_CascadeTools->a0xyError(moms[1],cascade_info->
getCovariance()[1],cascadeVertices[1],mainVertex);
251 for (
auto cascade_info : cascadeinfoContainer)
delete cascade_info;
252 for (
auto cascade_info_noConstr : cascadeinfoContainer_noConstr)
delete cascade_info_noConstr;
254 return StatusCode::SUCCESS;
258 m_vertexPsi1ContainerKey(
""),
259 m_vertexPsi2ContainerKey(
""),
260 m_cascadeOutputsKeys({
"PsiPlusPsiCascadeVtx1",
"PsiPlusPsiCascadeVtx2",
"PsiPlusPsiCascadeVtx3"}),
261 m_VxPrimaryCandidateName(
"PrimaryVertices"),
262 m_trackContainerName(
"InDetTrackParticles"),
263 m_eventInfo_key(
"EventInfo"),
264 m_jpsi1MassLower(0.0),
265 m_jpsi1MassUpper(20000.0),
266 m_jpsi2MassLower(0.0),
267 m_jpsi2MassUpper(20000.0),
268 m_diTrack1MassLower(-1.0),
269 m_diTrack1MassUpper(-1.0),
270 m_diTrack2MassLower(-1.0),
271 m_diTrack2MassUpper(-1.0),
272 m_psi1MassLower(0.0),
273 m_psi1MassUpper(25000.0),
274 m_psi2MassLower(0.0),
275 m_psi2MassUpper(25000.0),
277 m_MassUpper(31000.0),
279 m_vtx1Daug1MassHypo(-1),
280 m_vtx1Daug2MassHypo(-1),
281 m_vtx1Daug3MassHypo(-1),
282 m_vtx1Daug4MassHypo(-1),
284 m_vtx2Daug1MassHypo(-1),
285 m_vtx2Daug2MassHypo(-1),
286 m_vtx2Daug3MassHypo(-1),
287 m_vtx2Daug4MassHypo(-1),
296 m_constrJpsi1(
false),
297 m_constrJpsi2(
false),
298 m_constrDiTrk1(
false),
299 m_constrDiTrk2(
false),
300 m_chi2cut_Psi1(-1.0),
301 m_chi2cut_Psi2(-1.0),
303 m_removeDuplicatePairs(
false),
305 m_iVertexFitter(
"Trk::TrkVKalVrtFitter"),
306 m_pvRefitter(
"Analysis::PrimaryVertexRefitter",
this),
307 m_V0Tools(
"Trk::V0Tools"),
308 m_CascadeTools(
"DerivationFramework::CascadeTools")
316 declareProperty(
"RefPVContainerName", m_refPVContainerName =
"RefittedPrimaryVertices");
372 assert(cascadeinfoContainer!=
nullptr && cascadeinfoContainer_noConstr!=
nullptr);
378 std::vector<const xAOD::TrackParticle*> tracksJpsi1;
379 std::vector<const xAOD::TrackParticle*> tracksJpsi2;
380 std::vector<const xAOD::TrackParticle*> tracksDiTrk1;
381 std::vector<const xAOD::TrackParticle*> tracksDiTrk2;
382 std::vector<const xAOD::TrackParticle*> tracksPsi1;
383 std::vector<const xAOD::TrackParticle*> tracksPsi2;
384 std::vector<double> massesPsi1;
389 std::vector<double> massesPsi2;
404 std::vector<const xAOD::Vertex*> selectedPsi2Candidates;
405 for(
auto vxcItr=psi2Container.
cptr()->
cbegin(); vxcItr!=psi2Container.
cptr()->
cend(); ++vxcItr) {
418 double mass_psi2 =
m_V0Tools->invariantMass(*vxcItr, massesPsi2);
419 if (mass_psi2 < m_psi2MassLower || mass_psi2 >
m_psi2MassUpper)
continue;
422 TLorentzVector p4_mu1, p4_mu2;
423 p4_mu1.SetPtEtaPhiM( (*vxcItr)->trackParticle(0)->pt(),
424 (*vxcItr)->trackParticle(0)->eta(),
426 p4_mu2.SetPtEtaPhiM( (*vxcItr)->trackParticle(1)->pt(),
427 (*vxcItr)->trackParticle(1)->eta(),
429 double mass_jpsi2 = (p4_mu1 + p4_mu2).M();
430 if (mass_jpsi2 < m_jpsi2MassLower || mass_jpsi2 >
m_jpsi2MassUpper)
continue;
433 TLorentzVector p4_trk1, p4_trk2;
434 p4_trk1.SetPtEtaPhiM( (*vxcItr)->trackParticle(2)->pt(),
435 (*vxcItr)->trackParticle(2)->eta(),
437 p4_trk2.SetPtEtaPhiM( (*vxcItr)->trackParticle(3)->pt(),
438 (*vxcItr)->trackParticle(3)->eta(),
440 double mass_diTrk2 = (p4_trk1 + p4_trk2).M();
444 double chi2DOF = (*vxcItr)->chiSquared()/(*vxcItr)->numberDoF();
447 selectedPsi2Candidates.push_back(*vxcItr);
449 if(selectedPsi2Candidates.size()==0)
return StatusCode::SUCCESS;
452 std::vector<const xAOD::Vertex*> selectedPsi1Candidates;
453 for(
auto vxcItr=psi1Container.
cptr()->
cbegin(); vxcItr!=psi1Container.
cptr()->
cend(); ++vxcItr) {
466 double mass_psi1 =
m_V0Tools->invariantMass(*vxcItr,massesPsi1);
467 if(mass_psi1 < m_psi1MassLower || mass_psi1 >
m_psi1MassUpper)
continue;
470 TLorentzVector p4_mu1, p4_mu2;
471 p4_mu1.SetPtEtaPhiM( (*vxcItr)->trackParticle(0)->pt(),
472 (*vxcItr)->trackParticle(0)->eta(),
474 p4_mu2.SetPtEtaPhiM( (*vxcItr)->trackParticle(1)->pt(),
475 (*vxcItr)->trackParticle(1)->eta(),
477 double mass_jpsi1 = (p4_mu1 + p4_mu2).M();
478 if (mass_jpsi1 < m_jpsi1MassLower || mass_jpsi1 >
m_jpsi1MassUpper)
continue;
481 TLorentzVector p4_trk1, p4_trk2;
482 p4_trk1.SetPtEtaPhiM( (*vxcItr)->trackParticle(2)->pt(),
483 (*vxcItr)->trackParticle(2)->eta(),
485 p4_trk2.SetPtEtaPhiM( (*vxcItr)->trackParticle(3)->pt(),
486 (*vxcItr)->trackParticle(3)->eta(),
488 double mass_diTrk1 = (p4_trk1 + p4_trk2).M();
492 double chi2DOF = (*vxcItr)->chiSquared()/(*vxcItr)->numberDoF();
495 selectedPsi1Candidates.push_back(*vxcItr);
497 if(selectedPsi1Candidates.size()==0)
return StatusCode::SUCCESS;
499 std::vector<std::pair<const xAOD::Vertex*, const xAOD::Vertex*> > candidatePairs;
500 for(
auto psi1Itr=selectedPsi1Candidates.cbegin(); psi1Itr!=selectedPsi1Candidates.cend(); ++psi1Itr) {
502 for(
size_t i=0;
i<(*psi1Itr)->nTrackParticles();
i++) tracksPsi1.push_back((*psi1Itr)->trackParticle(
i));
503 for(
auto psi2Itr=selectedPsi2Candidates.cbegin(); psi2Itr!=selectedPsi2Candidates.cend(); ++psi2Itr) {
505 for(
size_t j=0; j<(*psi2Itr)->nTrackParticles(); j++) {
506 if(
std::find(tracksPsi1.cbegin(), tracksPsi1.cend(), (*psi2Itr)->trackParticle(j)) != tracksPsi1.cend()) {
skip =
true;
break; }
510 for(
size_t ic=0;
ic<candidatePairs.size();
ic++) {
513 if((psi1Vertex == *psi1Itr && psi2Vertex == *psi2Itr) || (psi1Vertex == *psi2Itr && psi2Vertex == *psi1Itr)) {
skip =
true;
break; }
517 candidatePairs.push_back(std::pair<const xAOD::Vertex*, const xAOD::Vertex*>(*psi1Itr,*psi2Itr));
521 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(); } );
523 candidatePairs.erase(candidatePairs.begin()+
m_maxCandidates, candidatePairs.end());
526 for(
size_t ic=0;
ic<candidatePairs.size();
ic++) {
532 if (tracksPsi1.size() != massesPsi1.size()) {
533 ATH_MSG_ERROR(
"Problems with Psi1 input: number of tracks or track mass inputs is not correct!");
537 if (tracksPsi2.size() != massesPsi2.size()) {
538 ATH_MSG_ERROR(
"Problems with Psi2 input: number of tracks or track mass inputs is not correct!");
544 tracksDiTrk1.clear();
552 tracksDiTrk2.clear();
558 TLorentzVector p4_moth;
577 std::vector<Trk::VertexID> vrtList;
586 vrtList.push_back(vID1);
594 vrtList.push_back(vID2);
596 std::vector<const xAOD::TrackParticle*>
tp;
tp.clear();
597 std::vector<double> tp_masses; tp_masses.clear();
600 std::vector<Trk::VertexID> cnstV; cnstV.clear();
606 std::vector<Trk::VertexID> cnstV; cnstV.clear();
612 std::vector<Trk::VertexID> cnstV; cnstV.clear();
618 std::vector<Trk::VertexID> cnstV; cnstV.clear();
628 for(
auto v :
result->vertices()) {
629 if(
v->nTrackParticles()==0) {
630 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
631 v->setTrackParticleLinks(nullLinkVector);
638 result->setSVOwnership(
true);
645 cascadeinfoContainer->push_back(
result.release());
655 std::vector<Trk::VertexID> vrtList_nc;
658 vrtList_nc.push_back(vID1_nc);
661 vrtList_nc.push_back(vID2_nc);
663 std::vector<const xAOD::TrackParticle*>
tp;
tp.clear();
664 std::vector<double> tp_masses; tp_masses.clear();
667 std::unique_ptr<Trk::VxCascadeInfo> result_nc(
m_iVertexFitter->fitCascade(*state));
669 if (result_nc !=
nullptr) {
671 if(
v->nTrackParticles()==0) {
672 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
673 v->setTrackParticleLinks(nullLinkVector);
681 cascadeinfoContainer_noConstr->push_back(result_nc.release());
683 else cascadeinfoContainer_noConstr->push_back(0);
685 else cascadeinfoContainer_noConstr->push_back(0);
689 return StatusCode::SUCCESS;