17 #include "HepPDT/ParticleDataTable.hh"
25 typedef std::vector<const xAOD::TrackParticle*>
TrackBag;
62 return StatusCode::SUCCESS;
68 std::vector<Trk::VxCascadeInfo*> cascadeinfoContainer;
69 constexpr
int topoN = 2;
70 std::array<xAOD::VertexContainer*, topoN> Vtxwritehandles;
71 std::array<xAOD::VertexAuxContainer*, topoN> Vtxwritehandlesaux;
74 for(
int i =0;
i<topoN;
i++){
77 Vtxwritehandles[
i]->setStore(Vtxwritehandlesaux[
i]);
90 if (pvContainer->
size()==0){
92 return StatusCode::RECOVERABLE;
94 primaryVertex = (*pvContainer)[0];
111 refPvContainer->setStore(refPvAuxContainer);
145 ATH_MSG_DEBUG(
"cascadeinfoContainer size " << cascadeinfoContainer.size());
157 return StatusCode::FAILURE;
178 const std::vector<xAOD::Vertex*> &cascadeVertices =
x->vertices();
179 if(cascadeVertices.size()!=topoN)
181 if(cascadeVertices[0] ==
nullptr || cascadeVertices[1] ==
nullptr)
ATH_MSG_ERROR(
"Error null vertex");
183 for(
int i =0;
i<topoN;
i++) Vtxwritehandles[
i]->push_back(cascadeVertices[
i]);
185 x->setSVOwnership(
false);
186 const auto mainVertex = cascadeVertices[1];
187 const std::vector< std::vector<TLorentzVector> > &moms =
x->getParticleMoms();
190 std::vector<const xAOD::Vertex*> verticestoLink;
191 verticestoLink.push_back(cascadeVertices[0]);
192 if(Vtxwritehandles[1] ==
nullptr)
ATH_MSG_ERROR(
"Vtxwritehandles[1] is null");
197 const xAOD::Vertex* jpsiVertex = BPhysPVCascadeTools::FindVertex<2>(jpsiContainer, cascadeVertices[1]);
198 ATH_MSG_DEBUG(
"1 pt Jpsi tracks " << cascadeVertices[1]->trackParticle(0)->
pt() <<
", " << cascadeVertices[1]->trackParticle(1)->
pt());
202 const xAOD::Vertex* dxVertex = BPhysPVCascadeTools::FindVertex<3>(dxContainer, cascadeVertices[0]);;
203 ATH_MSG_DEBUG(
"1 pt D_(s)+ tracks " << cascadeVertices[0]->trackParticle(0)->
pt() <<
", " << cascadeVertices[0]->trackParticle(1)->
pt() <<
", " << cascadeVertices[0]->trackParticle(2)->
pt());
207 std::vector<const xAOD::Vertex*> jpsiVerticestoLink;
208 if (jpsiVertex) jpsiVerticestoLink.push_back(jpsiVertex);
213 std::vector<const xAOD::Vertex*> dxVerticestoLink;
214 if (dxVertex) dxVerticestoLink.push_back(dxVertex);
226 std::vector<double> massesJpsi;
229 std::vector<double> massesDx;
238 std::vector<double> Masses;
260 PtErr_decor(*mainVertex) =
m_CascadeTools->pTError(moms[1],
x->getCovariance()[1]);
262 chi2_decor(*mainVertex) =
x->fitChi2();
263 ndof_decor(*mainVertex) =
x->nDoF();
267 TLorentzVector p4_mu1, p4_mu2;
274 massMumu = (p4_mu1 + p4_mu2).M();
276 MassMumu_decor(*mainVertex) = massMumu;
278 float massKX = 0., massKXpi = 0.;
280 TLorentzVector p4_h1, p4_h2, p4_h3;
299 massKX = (p4_h1 + p4_h2).M();
300 massKXpi = (p4_h1 + p4_h2 + p4_h3).M();
302 MassKX_svdecor(*mainVertex) = massKX;
303 MassKXpi_svdecor(*mainVertex) = massKXpi;
311 Mass_svdecor(*mainVertex) =
m_CascadeTools->invariantMass(moms[0]);
312 MassErr_svdecor(*mainVertex) =
m_CascadeTools->invariantMassError(moms[0],
x->getCovariance()[0]);
314 PtErr_svdecor(*mainVertex) =
m_CascadeTools->pTError(moms[0],
x->getCovariance()[0]);
315 Lxy_svdecor(*mainVertex) =
m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
316 LxyErr_svdecor(*mainVertex) =
m_CascadeTools->lxyError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
317 Tau_svdecor(*mainVertex) =
m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1]);
318 TauErr_svdecor(*mainVertex) =
m_CascadeTools->tauError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
322 <<
" chi2_1 " <<
m_V0Tools->chisq(cascadeVertices[0])
323 <<
" chi2_2 " <<
m_V0Tools->chisq(cascadeVertices[1])
327 <<
" error " <<
m_V0Tools->invariantMassError(cascadeVertices[0],massesDx)
328 <<
" mass_J " <<
m_V0Tools->invariantMass(cascadeVertices[1],massesJpsi)
329 <<
" error " <<
m_V0Tools->invariantMassError(cascadeVertices[1],massesJpsi));
333 double Mass_B_err =
m_CascadeTools->invariantMassError(moms[1],
x->getCovariance()[1]);
334 double Mass_D_err =
m_CascadeTools->invariantMassError(moms[0],
x->getCovariance()[0]);
336 ATH_MSG_DEBUG(
"Mass_B_err " << Mass_B_err <<
" Mass_D_err " << Mass_D_err);
338 double mprob_D =
m_CascadeTools->massProbability(mass_d,Mass_D,Mass_D_err);
339 ATH_MSG_DEBUG(
"mprob_B " << mprob_B <<
" mprob_D " << mprob_D);
342 <<
" Mass_d " <<
m_CascadeTools->invariantMass(moms[0],massesDx));
344 <<
" Mass_d_err " <<
m_CascadeTools->invariantMassError(moms[0],
x->getCovariance()[0],massesDx));
347 <<
" pt_dp " <<
m_V0Tools->pT(cascadeVertices[0]));
349 <<
" ptErr_d " <<
m_CascadeTools->pTError(moms[0],
x->getCovariance()[0])
350 <<
" ptErr_dp " <<
m_V0Tools->pTError(cascadeVertices[0]));
354 <<
" lxyErr_d " <<
m_CascadeTools->lxyError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
355 <<
" lxyErr_dp " <<
m_V0Tools->lxyError(cascadeVertices[0],cascadeVertices[1]));
357 <<
" tau_dp " <<
m_V0Tools->tau(cascadeVertices[0],cascadeVertices[1],massesDx));
359 <<
" tau_d " <<
m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1])
360 <<
" tau_D " <<
m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1],mass_d));
362 <<
" tauErr_d " <<
m_CascadeTools->tauError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
363 <<
" tauErr_dp " <<
m_V0Tools->tauError(cascadeVertices[0],cascadeVertices[1],massesDx));
365 <<
" TauErr_d " <<
m_CascadeTools->tauError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1],mass_d)
366 <<
" TauErr_dp " <<
m_V0Tools->tauError(cascadeVertices[0],cascadeVertices[1],massesDx,mass_d));
368 ATH_MSG_DEBUG(
"CascadeTools main vert wrt PV " <<
" CascadeTools SV " <<
" V0Tools SV");
370 <<
", " <<
m_CascadeTools->a0z(moms[0],cascadeVertices[0],cascadeVertices[1])
371 <<
", " <<
m_V0Tools->a0z(cascadeVertices[0],cascadeVertices[1]));
373 <<
", " <<
m_CascadeTools->a0zError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
374 <<
", " <<
m_V0Tools->a0zError(cascadeVertices[0],cascadeVertices[1]));
376 <<
", " <<
m_CascadeTools->a0xy(moms[0],cascadeVertices[0],cascadeVertices[1])
377 <<
", " <<
m_V0Tools->a0xy(cascadeVertices[0],cascadeVertices[1]));
379 <<
", " <<
m_CascadeTools->a0xyError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
380 <<
", " <<
m_V0Tools->a0xyError(cascadeVertices[0],cascadeVertices[1]));
382 <<
", " <<
m_CascadeTools->a0(moms[0],cascadeVertices[0],cascadeVertices[1])
383 <<
", " <<
m_V0Tools->a0(cascadeVertices[0],cascadeVertices[1]));
385 <<
", " <<
m_CascadeTools->a0Error(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
386 <<
", " <<
m_V0Tools->a0Error(cascadeVertices[0],cascadeVertices[1]));
389 ATH_MSG_DEBUG(
"X0 " << primaryVertex->
x() <<
" Y0 " << primaryVertex->
y() <<
" Z0 " << primaryVertex->
z());
392 ATH_MSG_DEBUG(
"Rxy0 wrt PV " <<
m_V0Tools->rxy(cascadeVertices[0],primaryVertex) <<
" RxyErr0 wrt PV " <<
m_V0Tools->rxyError(cascadeVertices[0],primaryVertex));
393 ATH_MSG_DEBUG(
"Rxy1 wrt PV " <<
m_V0Tools->rxy(cascadeVertices[1],primaryVertex) <<
" RxyErr1 wrt PV " <<
m_V0Tools->rxyError(cascadeVertices[1],primaryVertex));
394 ATH_MSG_DEBUG(
"number of covariance matrices " << (
x->getCovariance()).size());
399 for (
auto x : cascadeinfoContainer)
delete x;
401 return StatusCode::SUCCESS;
406 m_vertexContainerKey(
""),
407 m_vertexDxContainerKey(
""),
408 m_cascadeOutputsKeys{
"JpsiPlusDsCascadeVtx1",
"JpsiPlusDsCascadeVtx2" },
409 m_VxPrimaryCandidateName(
"PrimaryVertices"),
410 m_jpsiMassLower(0.0),
411 m_jpsiMassUpper(10000.0),
413 m_DxMassUpper(10000.0),
415 m_MassUpper(20000.0),
418 m_vtx0Daug1MassHypo(-1),
419 m_vtx0Daug2MassHypo(-1),
420 m_vtx1Daug1MassHypo(-1),
421 m_vtx1Daug2MassHypo(-1),
422 m_vtx1Daug3MassHypo(-1),
428 m_iVertexFitter(
"Trk::TrkVKalVrtFitter"),
429 m_pvRefitter(
"Analysis::PrimaryVertexRefitter"),
430 m_V0Tools(
"Trk::V0Tools"),
431 m_CascadeTools(
"DerivationFramework::CascadeTools")
433 declareProperty(
"JpsiVertices", m_vertexContainerKey);
434 declareProperty(
"DxVertices", m_vertexDxContainerKey);
435 declareProperty(
"VxPrimaryCandidateName", m_VxPrimaryCandidateName);
436 declareProperty(
"RefPVContainerName", m_refPVContainerName =
"RefittedPrimaryVertices");
437 declareProperty(
"JpsiMassLowerCut", m_jpsiMassLower);
438 declareProperty(
"JpsiMassUpperCut", m_jpsiMassUpper);
439 declareProperty(
"DxMassLowerCut", m_DxMassLower);
440 declareProperty(
"DxMassUpperCut", m_DxMassUpper);
441 declareProperty(
"MassLowerCut", m_MassLower);
442 declareProperty(
"MassUpperCut", m_MassUpper);
443 declareProperty(
"HypothesisName", m_hypoName =
"Bc");
444 declareProperty(
"Vtx0MassHypo", m_vtx0MassHypo);
445 declareProperty(
"Vtx1MassHypo", m_vtx1MassHypo);
446 declareProperty(
"Vtx0Daug1MassHypo", m_vtx0Daug1MassHypo);
447 declareProperty(
"Vtx0Daug2MassHypo", m_vtx0Daug2MassHypo);
448 declareProperty(
"Vtx1Daug1MassHypo", m_vtx1Daug1MassHypo);
449 declareProperty(
"Vtx1Daug2MassHypo", m_vtx1Daug2MassHypo);
450 declareProperty(
"Vtx1Daug3MassHypo", m_vtx1Daug3MassHypo);
451 declareProperty(
"JpsiMass", m_mass_jpsi);
452 declareProperty(
"DxHypothesis", m_Dx_pid);
453 declareProperty(
"ApplyDxMassConstraint", m_constrDx);
454 declareProperty(
"ApplyJpsiMassConstraint", m_constrJpsi);
455 declareProperty(
"Chi2Cut", m_chi2cut);
456 declareProperty(
"RefitPV", m_refitPV =
true);
457 declareProperty(
"MaxnPV", m_PV_max = 999);
458 declareProperty(
"MinNTracksInPV", m_PV_minNTracks = 0);
459 declareProperty(
"DoVertexType", m_DoVertexType = 7);
460 declareProperty(
"TrkVertexFitterTool", m_iVertexFitter);
461 declareProperty(
"PVRefitter", m_pvRefitter);
462 declareProperty(
"V0Tools", m_V0Tools);
463 declareProperty(
"CascadeTools", m_CascadeTools);
464 declareProperty(
"CascadeVertexCollections", m_cascadeOutputsKeys);
472 assert(cascadeinfoContainer!=
nullptr);
489 std::vector<const xAOD::TrackParticle*> tracksJpsi;
490 std::vector<const xAOD::TrackParticle*> tracksDx;
491 std::vector<const xAOD::TrackParticle*> tracksBc;
492 std::vector<double> massesJpsi;
495 std::vector<double> massesDx;
499 std::vector<double> massesDm;
503 std::vector<double> Masses;
509 std::vector<const xAOD::Vertex*> selectedJpsiCandidates;
510 for(
auto vxcItr=jpsiContainer->
cbegin(); vxcItr!=jpsiContainer->
cend(); ++vxcItr) {
516 if(!flagAcc1(*vtx))
continue;
520 double mass_Jpsi =
m_V0Tools->invariantMass(*vxcItr, massesJpsi);
522 ATH_MSG_DEBUG(
" Original Jpsi candidate rejected by the mass cut: mass = "
526 selectedJpsiCandidates.push_back(*vxcItr);
528 if(selectedJpsiCandidates.size()<1)
return StatusCode::SUCCESS;
531 std::vector<const xAOD::Vertex*> selectedDxCandidates;
532 for(
auto vxcItr=dxContainer->
cbegin(); vxcItr!=dxContainer->
cend(); ++vxcItr) {
539 if(!flagAcc1(*vtx))
continue;
549 if(!flagAcc1(*vtx)) isDp =
false;
552 if(!flagAcc2(*vtx)) isDm =
false;
554 if(!(isDp||isDm))
continue;
559 if(abs((*vxcItr)->trackParticle(0)->charge()+(*vxcItr)->trackParticle(1)->charge()+(*vxcItr)->trackParticle(2)->charge()) != 1){
560 ATH_MSG_DEBUG(
" Original D+ candidate rejected by the charge requirement: "
561 << (*vxcItr)->trackParticle(0)->charge() <<
", " << (*vxcItr)->trackParticle(1)->charge() <<
", " << (*vxcItr)->trackParticle(2)->charge() );
567 if(abs(
m_Dx_pid)==411 && (*vxcItr)->trackParticle(2)->charge()<0)
568 mass_D =
m_V0Tools->invariantMass(*vxcItr,massesDm);
570 mass_D =
m_V0Tools->invariantMass(*vxcItr,massesDx);
573 ATH_MSG_DEBUG(
" Original D_(s) candidate rejected by the mass cut: mass = "
580 TLorentzVector p4Kp_in, p4Km_in;
581 p4Kp_in.SetPtEtaPhiM( (*vxcItr)->trackParticle(0)->pt(),
582 (*vxcItr)->trackParticle(0)->eta(),
584 p4Km_in.SetPtEtaPhiM( (*vxcItr)->trackParticle(1)->pt(),
585 (*vxcItr)->trackParticle(1)->eta(),
587 double mass_phi = (p4Kp_in + p4Km_in).M();
589 if(mass_phi > 1200) {
590 ATH_MSG_DEBUG(
" Original phi candidate rejected by the mass cut: mass = " << mass_phi );
594 selectedDxCandidates.push_back(*vxcItr);
596 if(selectedDxCandidates.size()<1)
return StatusCode::SUCCESS;
600 for(
auto jpsiItr=selectedJpsiCandidates.cbegin(); jpsiItr!=selectedJpsiCandidates.cend(); ++jpsiItr) {
602 size_t jpsiTrkNum = (*jpsiItr)->nTrackParticles();
604 for(
unsigned int it=0;
it<jpsiTrkNum;
it++) tracksJpsi.push_back((*jpsiItr)->trackParticle(
it));
606 if (tracksJpsi.size() != 2 || massesJpsi.size() != 2 ) {
611 for(
auto dxItr=selectedDxCandidates.cbegin(); dxItr!=selectedDxCandidates.cend(); ++dxItr) {
614 if(
std::find(tracksJpsi.cbegin(), tracksJpsi.cend(), (*dxItr)->trackParticle(0)) != tracksJpsi.cend())
continue;
615 if(
std::find(tracksJpsi.cbegin(), tracksJpsi.cend(), (*dxItr)->trackParticle(1)) != tracksJpsi.cend())
continue;
616 if(
std::find(tracksJpsi.cbegin(), tracksJpsi.cend(), (*dxItr)->trackParticle(2)) != tracksJpsi.cend())
continue;
618 size_t dxTrkNum = (*dxItr)->nTrackParticles();
620 for(
unsigned int it=0;
it<dxTrkNum;
it++) tracksDx.push_back((*dxItr)->trackParticle(
it));
621 if (tracksDx.size() != 3 || massesDx.size() != 3 ) {
625 ATH_MSG_DEBUG(
"using tracks" << tracksJpsi[0] <<
", " << tracksJpsi[1] <<
", " << tracksDx[0] <<
", " << tracksDx[1] <<
", " << tracksDx[2]);
627 for(
unsigned int it=0;
it<jpsiTrkNum;
it++) tracksBc.push_back((*jpsiItr)->trackParticle(
it));
628 for(
unsigned int it=0;
it<dxTrkNum;
it++) tracksBc.push_back((*dxItr)->trackParticle(
it));
638 std::vector<Trk::VertexID> vrtList;
642 if(abs(
m_Dx_pid)==411 && (*dxItr)->trackParticle(2)->charge()<0)
647 if(abs(
m_Dx_pid)==411 && (*dxItr)->trackParticle(2)->charge()<0)
652 vrtList.push_back(vID);
656 std::vector<Trk::VertexID> cnstV;
670 << ((
result->vertices())[0])->trackParticle(1) <<
", "
671 << ((
result->vertices())[0])->trackParticle(2) <<
", "
672 << ((
result->vertices())[1])->trackParticle(0) <<
", "
673 << ((
result->vertices())[1])->trackParticle(1));
675 result->setSVOwnership(
true);
682 const std::vector< std::vector<TLorentzVector> > &moms =
result->getParticleMoms();
686 cascadeinfoContainer->push_back(
result.release());
698 ATH_MSG_DEBUG(
"cascadeinfoContainer size " << cascadeinfoContainer->size());
700 return StatusCode::SUCCESS;