11 #include "GaudiKernel/IPartPropSvc.h"
18 #include "HepPDT/ParticleDataTable.hh"
26 typedef std::vector<const xAOD::TrackParticle*>
TrackBag;
43 IPartPropSvc* partPropSvc =
nullptr;
44 ATH_CHECK( service(
"PartPropSvc", partPropSvc,
true) );
45 auto pdt = partPropSvc->PDT();
58 return StatusCode::SUCCESS;
64 std::vector<Trk::VxCascadeInfo*> cascadeinfoContainer;
65 constexpr
int topoN = 2;
66 std::array<xAOD::VertexContainer*, topoN> Vtxwritehandles;
67 std::array<xAOD::VertexAuxContainer*, topoN> Vtxwritehandlesaux;
70 for(
int i =0;
i<topoN;
i++){
73 Vtxwritehandles[
i]->setStore(Vtxwritehandlesaux[
i]);
86 if (pvContainer->
size()==0){
88 return StatusCode::RECOVERABLE;
90 primaryVertex = (*pvContainer)[0];
107 refPvContainer->setStore(refPvAuxContainer);
142 ATH_MSG_DEBUG(
"cascadeinfoContainer size " << cascadeinfoContainer.size());
154 return StatusCode::FAILURE;
175 const std::vector<xAOD::Vertex*> &cascadeVertices =
x->vertices();
176 if(cascadeVertices.size()!=topoN)
178 if(cascadeVertices[0] ==
nullptr || cascadeVertices[1] ==
nullptr)
ATH_MSG_ERROR(
"Error null vertex");
180 for(
int i =0;
i<topoN;
i++) Vtxwritehandles[
i]->push_back(cascadeVertices[
i]);
182 x->setSVOwnership(
false);
183 const auto mainVertex = cascadeVertices[1];
184 const std::vector< std::vector<TLorentzVector> > &moms =
x->getParticleMoms();
187 std::vector<const xAOD::Vertex*> verticestoLink;
188 verticestoLink.push_back(cascadeVertices[0]);
189 if(Vtxwritehandles[1] ==
nullptr)
ATH_MSG_ERROR(
"Vtxwritehandles[1] is null");
194 const xAOD::Vertex* jpsipiVertex = BPhysPVCascadeTools::FindVertex<3>(jpsipiContainer, cascadeVertices[1]);
195 ATH_MSG_DEBUG(
"1 pt Jpsi+pi tracks " << cascadeVertices[1]->trackParticle(0)->
pt() <<
", " << cascadeVertices[1]->trackParticle(1)->
pt() <<
", " << cascadeVertices[1]->trackParticle(2)->
pt());
199 const xAOD::Vertex* d0Vertex = BPhysPVCascadeTools::FindVertex<2>(d0Container, cascadeVertices[0]);;
200 ATH_MSG_DEBUG(
"1 pt D0 tracks " << cascadeVertices[0]->trackParticle(0)->
pt() <<
", " << cascadeVertices[0]->trackParticle(1)->
pt());
204 std::vector<const xAOD::Vertex*> jpsipiVerticestoLink;
205 if (jpsipiVertex) jpsipiVerticestoLink.push_back(jpsipiVertex);
210 std::vector<const xAOD::Vertex*> d0VerticestoLink;
211 if (d0Vertex) d0VerticestoLink.push_back(d0Vertex);
223 std::vector<double> massesJpsipi;
227 std::vector<double> massesD0;
235 std::vector<double> Masses;
257 PtErr_decor(*mainVertex) =
m_CascadeTools->pTError(moms[1],
x->getCovariance()[1]);
259 chi2_decor(*mainVertex) =
x->fitChi2();
260 ndof_decor(*mainVertex) =
x->nDoF();
264 TLorentzVector p4_mu1, p4_mu2;
271 massMumu = (p4_mu1 + p4_mu2).M();
273 MassMumu_decor(*mainVertex) = massMumu;
277 TLorentzVector p4_ka, p4_pi;
293 massKpi = (p4_ka + p4_pi).M();
295 MassKpi_svdecor(*mainVertex) = massKpi;
296 MassJpsi_decor(*mainVertex) = (moms[1][0] + moms[1][1]).M();
297 MassPiD0_decor(*mainVertex) = (moms[1][2] + moms[1][3]).M();
305 Mass_svdecor(*mainVertex) =
m_CascadeTools->invariantMass(moms[0]);
306 MassErr_svdecor(*mainVertex) =
m_CascadeTools->invariantMassError(moms[0],
x->getCovariance()[0]);
308 PtErr_svdecor(*mainVertex) =
m_CascadeTools->pTError(moms[0],
x->getCovariance()[0]);
309 Lxy_svdecor(*mainVertex) =
m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
310 LxyErr_svdecor(*mainVertex) =
m_CascadeTools->lxyError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
311 Tau_svdecor(*mainVertex) =
m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1]);
312 TauErr_svdecor(*mainVertex) =
m_CascadeTools->tauError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
316 <<
" chi2_1 " <<
m_V0Tools->chisq(cascadeVertices[0])
317 <<
" chi2_2 " <<
m_V0Tools->chisq(cascadeVertices[1])
321 <<
" error " <<
m_V0Tools->invariantMassError(cascadeVertices[0],massesD0)
322 <<
" mass_J " <<
m_V0Tools->invariantMass(cascadeVertices[1],massesJpsipi)
323 <<
" error " <<
m_V0Tools->invariantMassError(cascadeVertices[1],massesJpsipi));
327 double Mass_B_err =
m_CascadeTools->invariantMassError(moms[1],
x->getCovariance()[1]);
328 double Mass_D0_err =
m_CascadeTools->invariantMassError(moms[0],
x->getCovariance()[0]);
329 ATH_MSG_DEBUG(
"Mass_B " << Mass_B <<
" Mass_D0 " << Mass_D0);
330 ATH_MSG_DEBUG(
"Mass_B_err " << Mass_B_err <<
" Mass_D0_err " << Mass_D0_err);
332 double mprob_D0 =
m_CascadeTools->massProbability(mass_d0,Mass_D0,Mass_D0_err);
333 ATH_MSG_DEBUG(
"mprob_B " << mprob_B <<
" mprob_D0 " << mprob_D0);
336 <<
" Mass_d0 " <<
m_CascadeTools->invariantMass(moms[0],massesD0));
338 <<
" Mass_d0_err " <<
m_CascadeTools->invariantMassError(moms[0],
x->getCovariance()[0],massesD0));
341 <<
" pt_d0 " <<
m_V0Tools->pT(cascadeVertices[0]));
343 <<
" ptErr_d " <<
m_CascadeTools->pTError(moms[0],
x->getCovariance()[0])
344 <<
" ptErr_d0 " <<
m_V0Tools->pTError(cascadeVertices[0]));
348 <<
" lxyErr_d " <<
m_CascadeTools->lxyError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
349 <<
" lxyErr_d0 " <<
m_V0Tools->lxyError(cascadeVertices[0],cascadeVertices[1]));
351 <<
" tau_d0 " <<
m_V0Tools->tau(cascadeVertices[0],cascadeVertices[1],massesD0));
353 <<
" tau_d " <<
m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1])
354 <<
" tau_D " <<
m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1],mass_d0));
356 <<
" tauErr_d " <<
m_CascadeTools->tauError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
357 <<
" tauErr_d0 " <<
m_V0Tools->tauError(cascadeVertices[0],cascadeVertices[1],massesD0));
359 <<
" TauErr_d " <<
m_CascadeTools->tauError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1],mass_d0)
360 <<
" TauErr_d0 " <<
m_V0Tools->tauError(cascadeVertices[0],cascadeVertices[1],massesD0,mass_d0));
362 ATH_MSG_DEBUG(
"CascadeTools main vert wrt PV " <<
" CascadeTools SV " <<
" V0Tools SV");
364 <<
", " <<
m_CascadeTools->a0z(moms[0],cascadeVertices[0],cascadeVertices[1])
365 <<
", " <<
m_V0Tools->a0z(cascadeVertices[0],cascadeVertices[1]));
367 <<
", " <<
m_CascadeTools->a0zError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
368 <<
", " <<
m_V0Tools->a0zError(cascadeVertices[0],cascadeVertices[1]));
370 <<
", " <<
m_CascadeTools->a0xy(moms[0],cascadeVertices[0],cascadeVertices[1])
371 <<
", " <<
m_V0Tools->a0xy(cascadeVertices[0],cascadeVertices[1]));
373 <<
", " <<
m_CascadeTools->a0xyError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
374 <<
", " <<
m_V0Tools->a0xyError(cascadeVertices[0],cascadeVertices[1]));
376 <<
", " <<
m_CascadeTools->a0(moms[0],cascadeVertices[0],cascadeVertices[1])
377 <<
", " <<
m_V0Tools->a0(cascadeVertices[0],cascadeVertices[1]));
379 <<
", " <<
m_CascadeTools->a0Error(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
380 <<
", " <<
m_V0Tools->a0Error(cascadeVertices[0],cascadeVertices[1]));
383 ATH_MSG_DEBUG(
"X0 " << primaryVertex->
x() <<
" Y0 " << primaryVertex->
y() <<
" Z0 " << primaryVertex->
z());
386 ATH_MSG_DEBUG(
"Rxy0 wrt PV " <<
m_V0Tools->rxy(cascadeVertices[0],primaryVertex) <<
" RxyErr0 wrt PV " <<
m_V0Tools->rxyError(cascadeVertices[0],primaryVertex));
387 ATH_MSG_DEBUG(
"Rxy1 wrt PV " <<
m_V0Tools->rxy(cascadeVertices[1],primaryVertex) <<
" RxyErr1 wrt PV " <<
m_V0Tools->rxyError(cascadeVertices[1],primaryVertex));
388 ATH_MSG_DEBUG(
"number of covariance matrices " << (
x->getCovariance()).size());
393 for (
auto x : cascadeinfoContainer)
delete x;
395 return StatusCode::SUCCESS;
400 m_vertexContainerKey(
""),
401 m_vertexD0ContainerKey(
""),
402 m_cascadeOutputsKeys{
"JpsiPlusDpstCascadeVtx1",
"JpsiPlusDpstCascadeVtx2" },
403 m_VxPrimaryCandidateName(
"PrimaryVertices"),
404 m_jpsiMassLower(0.0),
405 m_jpsiMassUpper(10000.0),
406 m_jpsipiMassLower(0.0),
407 m_jpsipiMassUpper(10000.0),
409 m_D0MassUpper(10000.0),
411 m_DstMassUpper(10000.0),
413 m_MassUpper(20000.0),
416 m_vtx0Daug1MassHypo(-1),
417 m_vtx0Daug2MassHypo(-1),
418 m_vtx0Daug3MassHypo(-1),
419 m_vtx1Daug1MassHypo(-1),
420 m_vtx1Daug2MassHypo(-1),
426 m_iVertexFitter(
"Trk::TrkVKalVrtFitter"),
427 m_pvRefitter(
"Analysis::PrimaryVertexRefitter"),
428 m_V0Tools(
"Trk::V0Tools"),
429 m_CascadeTools(
"DerivationFramework::CascadeTools")
434 declareProperty(
"RefPVContainerName", m_refPVContainerName =
"RefittedPrimaryVertices");
474 assert(cascadeinfoContainer!=
nullptr);
489 std::vector<const xAOD::TrackParticle*> tracksJpsipi;
490 std::vector<const xAOD::TrackParticle*> tracksJpsi;
491 std::vector<const xAOD::TrackParticle*> tracksD0;
492 std::vector<const xAOD::TrackParticle*> tracksBc;
493 std::vector<double> massesJpsipi;
497 std::vector<double> massesD0;
500 std::vector<double> massesD0b;
503 std::vector<double> Masses;
510 std::vector<const xAOD::Vertex*> selectedJpsipiCandidates;
511 for(
auto vxcItr=jpsipiContainer->
cbegin(); vxcItr!=jpsipiContainer->
cend(); ++vxcItr) {
517 if(!flagAcc1(*vtx))
continue;
521 TLorentzVector p4Mup_in, p4Mum_in;
522 p4Mup_in.SetPtEtaPhiM((*vxcItr)->trackParticle(0)->pt(),
523 (*vxcItr)->trackParticle(0)->eta(),
525 p4Mum_in.SetPtEtaPhiM((*vxcItr)->trackParticle(1)->pt(),
526 (*vxcItr)->trackParticle(1)->eta(),
528 double mass_Jpsi = (p4Mup_in + p4Mum_in).M();
531 ATH_MSG_DEBUG(
" Original Jpsi candidate rejected by the mass cut: mass = "
537 double mass_Jpsipi =
m_V0Tools->invariantMass(*vxcItr, massesJpsipi);
540 ATH_MSG_DEBUG(
" Original Jpsipi candidate rejected by the mass cut: mass = "
545 selectedJpsipiCandidates.push_back(*vxcItr);
547 if(selectedJpsipiCandidates.size()<1)
return StatusCode::SUCCESS;
550 std::vector<const xAOD::Vertex*> selectedD0Candidates;
551 for(
auto vxcItr=d0Container->
cbegin(); vxcItr!=d0Container->
cend(); ++vxcItr) {
560 if(!flagAcc1(*vtx)) isD0 =
false;
563 if(!flagAcc2(*vtx)) isD0b =
false;
565 if(!(isD0||isD0b))
continue;
568 if ((*vxcItr)->trackParticle(0)->charge() != 1 || (*vxcItr)->trackParticle(1)->charge() != -1) {
569 ATH_MSG_DEBUG(
" Original D0/D0-bar candidate rejected by the charge requirement: "
570 << (*vxcItr)->trackParticle(0)->charge() <<
", " << (*vxcItr)->trackParticle(1)->charge() );
575 double mass_D0 =
m_V0Tools->invariantMass(*vxcItr,massesD0);
576 double mass_D0b =
m_V0Tools->invariantMass(*vxcItr,massesD0b);
577 ATH_MSG_DEBUG(
"D0 mass " << mass_D0 <<
", D0b mass "<<mass_D0b);
579 ATH_MSG_DEBUG(
" Original D0 candidate rejected by the mass cut: mass = "
585 selectedD0Candidates.push_back(*vxcItr);
587 if(selectedD0Candidates.size()<1)
return StatusCode::SUCCESS;
591 for(
auto jpsipiItr=selectedJpsipiCandidates.cbegin(); jpsipiItr!=selectedJpsipiCandidates.cend(); ++jpsipiItr) {
593 size_t jpsipiTrkNum = (*jpsipiItr)->nTrackParticles();
594 tracksJpsipi.clear();
596 for(
unsigned int it=0;
it<jpsipiTrkNum;
it++) tracksJpsipi.push_back((*jpsipiItr)->trackParticle(
it));
597 for(
unsigned int it=0;
it<jpsipiTrkNum-1;
it++) tracksJpsi.push_back((*jpsipiItr)->trackParticle(
it));
599 if (tracksJpsipi.size() != 3 || massesJpsipi.size() != 3 ) {
604 if(abs(
m_Dx_pid)==421 && (*jpsipiItr)->trackParticle(2)->charge()==-1) tagD0 =
false;
606 TLorentzVector p4_pi1;
607 p4_pi1.SetPtEtaPhiM((*jpsipiItr)->trackParticle(2)->pt(),
608 (*jpsipiItr)->trackParticle(2)->eta(),
612 for(
auto d0Itr=selectedD0Candidates.cbegin(); d0Itr!=selectedD0Candidates.cend(); ++d0Itr) {
615 if(
std::find(tracksJpsipi.cbegin(), tracksJpsipi.cend(), (*d0Itr)->trackParticle(0)) != tracksJpsipi.cend())
continue;
616 if(
std::find(tracksJpsipi.cbegin(), tracksJpsipi.cend(), (*d0Itr)->trackParticle(1)) != tracksJpsipi.cend())
continue;
619 TLorentzVector p4_ka, p4_pi2;
621 p4_pi2.SetPtEtaPhiM((*d0Itr)->trackParticle(0)->pt(),
622 (*d0Itr)->trackParticle(0)->eta(),
624 p4_ka.SetPtEtaPhiM( (*d0Itr)->trackParticle(1)->pt(),
625 (*d0Itr)->trackParticle(1)->eta(),
628 p4_pi2.SetPtEtaPhiM((*d0Itr)->trackParticle(1)->pt(),
629 (*d0Itr)->trackParticle(1)->eta(),
631 p4_ka.SetPtEtaPhiM( (*d0Itr)->trackParticle(0)->pt(),
632 (*d0Itr)->trackParticle(0)->eta(),
636 double mass_Dst= (p4_pi1 + p4_ka + p4_pi2).M();
639 ATH_MSG_DEBUG(
" Original D*+/- candidate rejected by the mass cut: mass = "
644 size_t d0TrkNum = (*d0Itr)->nTrackParticles();
646 for(
unsigned int it=0;
it<d0TrkNum;
it++) tracksD0.push_back((*d0Itr)->trackParticle(
it));
647 if (tracksD0.size() != 2 || massesD0.size() != 2 ) {
651 ATH_MSG_DEBUG(
"using tracks" << tracksJpsipi[0] <<
", " << tracksJpsipi[1] <<
", " << tracksJpsipi[2] <<
", " << tracksD0[0] <<
", " << tracksD0[1]);
652 ATH_MSG_DEBUG(
"Charge of Jpsi+pi tracks: "<<(*jpsipiItr)->trackParticle(0)->charge()<<
", "<<(*jpsipiItr)->trackParticle(1)->charge()<<
", "<<(*jpsipiItr)->trackParticle(2)->charge());
653 ATH_MSG_DEBUG(
"Charge of D0 tracks: "<<(*d0Itr)->trackParticle(0)->charge()<<
", "<<(*d0Itr)->trackParticle(1)->charge());
656 for(
unsigned int it=0;
it<jpsipiTrkNum;
it++) tracksBc.push_back((*jpsipiItr)->trackParticle(
it));
657 for(
unsigned int it=0;
it<d0TrkNum;
it++) tracksBc.push_back((*d0Itr)->trackParticle(
it));
668 std::vector<Trk::VertexID> vrtList;
672 if(tagD0) vID =
m_iVertexFitter->startVertex(tracksD0,massesD0,*state,mass_d0);
673 else vID =
m_iVertexFitter->startVertex(tracksD0,massesD0b,*state,mass_d0);
675 if(tagD0) vID =
m_iVertexFitter->startVertex(tracksD0,massesD0,*state);
678 vrtList.push_back(vID);
682 std::vector<Trk::VertexID> cnstV;
697 << ((
result->vertices())[0])->trackParticle(1) <<
", "
698 << ((
result->vertices())[1])->trackParticle(0) <<
", "
699 << ((
result->vertices())[1])->trackParticle(1) <<
", "
700 << ((
result->vertices())[1])->trackParticle(2));
702 result->setSVOwnership(
true);
709 const std::vector< std::vector<TLorentzVector> > &moms =
result->getParticleMoms();
713 cascadeinfoContainer->push_back(
result.release());
725 ATH_MSG_DEBUG(
"cascadeinfoContainer size " << cascadeinfoContainer->size());
727 return StatusCode::SUCCESS;