17 #include "HepPDT/ParticleDataTable.hh"
33 typedef std::vector<const xAOD::TrackParticle*>
TrackBag;
49 return StatusCode::FAILURE;
91 return StatusCode::SUCCESS;
97 std::vector<Trk::VxCascadeInfo*> cascadeinfoContainer;
98 constexpr
int topoN = 2;
99 std::array<xAOD::VertexContainer*, topoN> Vtxwritehandles;
100 std::array<xAOD::VertexAuxContainer*, topoN> Vtxwritehandlesaux;
103 for(
int i =0;
i<topoN;
i++){
106 Vtxwritehandles[
i]->setStore(Vtxwritehandlesaux[
i]);
119 if (pvContainer->
size()==0){
121 return StatusCode::RECOVERABLE;
123 primaryVertex = (*pvContainer)[0];
140 refPvContainer->setStore(refPvAuxContainer);
150 return StatusCode::FAILURE;
198 ATH_MSG_DEBUG(
"cascadeinfoContainer size " << cascadeinfoContainer.size());
233 const std::vector<xAOD::Vertex*> &cascadeVertices =
x->vertices();
235 if(cascadeVertices.size()!=topoN)
ATH_MSG_ERROR(
"Incorrect number of vertices");
236 if(cascadeVertices[0] ==
nullptr || cascadeVertices[1] ==
nullptr)
ATH_MSG_ERROR(
"Error null vertex");
238 for(
int i =0;
i<topoN;
i++) Vtxwritehandles[
i]->push_back(cascadeVertices[
i]);
240 x->setSVOwnership(
false);
241 const auto mainVertex = cascadeVertices[1];
243 const std::vector< std::vector<TLorentzVector> > &moms =
x->getParticleMoms();
247 std::vector<const xAOD::Vertex*> verticestoLink;
248 verticestoLink.push_back(cascadeVertices[0]);
249 if(Vtxwritehandles[1] ==
nullptr)
ATH_MSG_ERROR(
"Vtxwritehandles[1] is null");
255 typedef std::vector<const xAOD::Muon*>
MuonBag;
256 MuonBag selectedMuons; selectedMuons.clear();
260 if (muonTrk == cascadeVertices[1]->trackParticle(0)) selectedMuons.push_back(
mu);
270 auto muItr = selectedMuons.begin();
271 for(; muItr != selectedMuons.end(); muItr++) {
278 if( !muLink.
isValid() )
continue;
281 preMuLinks.push_back( muLink );
286 MuonsLinksDecor(*cascadeVertices[1]) = preMuLinks;
292 const xAOD::Vertex* dxVertex = BPhysPVCascadeTools::FindVertex<3>(dxContainer, cascadeVertices[0]);
293 ATH_MSG_DEBUG(
"1 pt D_(s)+/Lambda_c+ tracks " << cascadeVertices[1]->trackParticle(0)->
pt()<<
", "
294 << cascadeVertices[0]->trackParticle(0)->
pt()<<
", "<< cascadeVertices[0]->trackParticle(1)->
pt()<<
", "
295 << cascadeVertices[0]->trackParticle(2)->
pt());
299 std::vector<const xAOD::Vertex*> dxVerticestoLink;
300 if (dxVertex) dxVerticestoLink.push_back(dxVertex);
303 ATH_MSG_ERROR(
"Error decorating with D_(s)+/Lambda_c+ vertices");
314 std::vector<double> massesMu;
316 std::vector<double> massesDx;
325 std::vector<double> Masses;
346 PtErr_decor(*mainVertex) =
m_CascadeTools->pTError(moms[1],
x->getCovariance()[1]);
348 chi2_decor(*mainVertex) =
x->fitChi2();
349 ndof_decor(*mainVertex) =
x->nDoF();
350 Dx_chi2_svdecor(*mainVertex) = cascadeVertices[0]->chiSquared();
351 Dx_ndof_svdecor(*mainVertex) = cascadeVertices[0]->numberDoF();
354 float muM = 0., mupt = 0., mueta = 0.;
355 if (selectedMuons.size()==1) {
356 TLorentzVector p4_mu1;
357 p4_mu1.SetPtEtaPhiM(selectedMuons.at(0)->pt(),
358 selectedMuons.at(0)->eta(),
362 mueta = p4_mu1.Eta();
364 MuMass_decor(*mainVertex) = muM;
365 MuPt_decor(*mainVertex) = mupt;
366 MuEta_decor(*mainVertex) = mueta;
367 ChargeMu_decor(*mainVertex) = selectedMuons.at(0)->charge();
368 MuChi2_decor(*mainVertex) = cascadeVertices[1]->trackParticle(0)->chiSquared();
369 MunDoF_decor(*mainVertex) = cascadeVertices[1]->trackParticle(0)->numberDoF();
372 std::vector< Trk::VxTrackAtVertex > trkAtB = cascadeVertices[1]->vxTrackAtVertex();
373 MuChi2B_decor(*mainVertex) = trkAtB.at(0).trackQuality().chiSquared();
374 MunDoFB_decor(*mainVertex) = trkAtB.at(0).trackQuality().numberDoF();
379 float massKX1X2 = 0.;
380 float RapidityKX1X2 = 0.;
382 TLorentzVector p4_h1, p4_h2, p4_h3;
408 massKX1 = (p4_h1 + p4_h2).M();
409 massKX1X2 = (p4_h1 + p4_h2 + p4_h3).M();
410 RapidityKX1X2 = (p4_h1 + p4_h2 + p4_h3).Rapidity();
412 massKX1_svdecor(*mainVertex) = massKX1;
413 massKX1X2_svdecor(*mainVertex) = massKX1X2;
414 RapidityKX1X2_svdecor(*mainVertex) = RapidityKX1X2;
423 Mass_svdecor(*mainVertex) =
m_CascadeTools->invariantMass(moms[0]);
424 MassErr_svdecor(*mainVertex) =
m_CascadeTools->invariantMassError(moms[0],
x->getCovariance()[0]);
426 PtErr_svdecor(*mainVertex) =
m_CascadeTools->pTError(moms[0],
x->getCovariance()[0]);
427 Lxy_svdecor(*mainVertex) =
m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
428 LxyErr_svdecor(*mainVertex) =
m_CascadeTools->lxyError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
429 Tau_svdecor(*mainVertex) =
m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1]);
430 TauErr_svdecor(*mainVertex) =
m_CascadeTools->tauError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
434 <<
" chi2_1 " <<
m_V0Tools->chisq(cascadeVertices[0])
435 <<
" chi2_2 " <<
m_V0Tools->chisq(cascadeVertices[1])
439 <<
" error " <<
m_V0Tools->invariantMassError(cascadeVertices[0],massesDx)
440 <<
" mass_J " <<
m_V0Tools->invariantMass(cascadeVertices[1],massesMu)
441 <<
" error " <<
m_V0Tools->invariantMassError(cascadeVertices[1],massesMu));
445 double Mass_B_err =
m_CascadeTools->invariantMassError(moms[1],
x->getCovariance()[1]);
446 double Mass_D_err =
m_CascadeTools->invariantMassError(moms[0],
x->getCovariance()[0]);
448 ATH_MSG_DEBUG(
"Mass_B_err " << Mass_B_err <<
" Mass_D_err " << Mass_D_err);
450 double mprob_D =
m_CascadeTools->massProbability(mass_d,Mass_D,Mass_D_err);
451 ATH_MSG_DEBUG(
"mprob_B " << mprob_B <<
" mprob_D " << mprob_D);
454 <<
" Mass_d " <<
m_CascadeTools->invariantMass(moms[0],massesDx));
456 <<
" Mass_d_err " <<
m_CascadeTools->invariantMassError(moms[0],
x->getCovariance()[0],massesDx));
459 <<
" pt_dp " <<
m_V0Tools->pT(cascadeVertices[0]));
461 <<
" ptErr_d " <<
m_CascadeTools->pTError(moms[0],
x->getCovariance()[0])
462 <<
" ptErr_dp " <<
m_V0Tools->pTError(cascadeVertices[0]));
466 <<
" lxyErr_d " <<
m_CascadeTools->lxyError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
467 <<
" lxyErr_dp " <<
m_V0Tools->lxyError(cascadeVertices[0],cascadeVertices[1]));
469 <<
" tau_dp " <<
m_V0Tools->tau(cascadeVertices[0],cascadeVertices[1],massesDx));
471 <<
" tau_d " <<
m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1])
472 <<
" tau_D " <<
m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1],mass_d));
474 <<
" tauErr_d " <<
m_CascadeTools->tauError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
475 <<
" tauErr_dp " <<
m_V0Tools->tauError(cascadeVertices[0],cascadeVertices[1],massesDx));
477 <<
" TauErr_d " <<
m_CascadeTools->tauError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1],mass_d)
478 <<
" TauErr_dp " <<
m_V0Tools->tauError(cascadeVertices[0],cascadeVertices[1],massesDx,mass_d));
480 ATH_MSG_DEBUG(
"CascadeTools main vert wrt PV " <<
" CascadeTools SV " <<
" V0Tools SV");
482 <<
", " <<
m_CascadeTools->a0z(moms[0],cascadeVertices[0],cascadeVertices[1])
483 <<
", " <<
m_V0Tools->a0z(cascadeVertices[0],cascadeVertices[1]));
485 <<
", " <<
m_CascadeTools->a0zError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
486 <<
", " <<
m_V0Tools->a0zError(cascadeVertices[0],cascadeVertices[1]));
488 <<
", " <<
m_CascadeTools->a0xy(moms[0],cascadeVertices[0],cascadeVertices[1])
489 <<
", " <<
m_V0Tools->a0xy(cascadeVertices[0],cascadeVertices[1]));
491 <<
", " <<
m_CascadeTools->a0xyError(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
492 <<
", " <<
m_V0Tools->a0xyError(cascadeVertices[0],cascadeVertices[1]));
494 <<
", " <<
m_CascadeTools->a0(moms[0],cascadeVertices[0],cascadeVertices[1])
495 <<
", " <<
m_V0Tools->a0(cascadeVertices[0],cascadeVertices[1]));
497 <<
", " <<
m_CascadeTools->a0Error(moms[0],
x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
498 <<
", " <<
m_V0Tools->a0Error(cascadeVertices[0],cascadeVertices[1]));
501 ATH_MSG_DEBUG(
"X0 " << primaryVertex->
x() <<
" Y0 " << primaryVertex->
y() <<
" Z0 " << primaryVertex->
z());
504 ATH_MSG_DEBUG(
"Rxy0 wrt PV " <<
m_V0Tools->rxy(cascadeVertices[0],primaryVertex) <<
" RxyErr0 wrt PV " <<
m_V0Tools->rxyError(cascadeVertices[0],primaryVertex));
505 ATH_MSG_DEBUG(
"Rxy1 wrt PV " <<
m_V0Tools->rxy(cascadeVertices[1],primaryVertex) <<
" RxyErr1 wrt PV " <<
m_V0Tools->rxyError(cascadeVertices[1],primaryVertex));
506 ATH_MSG_DEBUG(
"number of covariance matrices " << (
x->getCovariance()).size());
512 for (
auto x : cascadeinfoContainer)
delete x;
514 return StatusCode::SUCCESS;
519 m_vertexDxContainerKey(
""),
520 m_cascadeOutputsKeys{
"MuPlusDsCascadeVtx1",
"MuPlusDsCascadeVtx2" },
521 m_VxPrimaryCandidateName(
"PrimaryVertices"),
523 m_DxMassUpper(10000.0),
525 m_MassUpper(20000.0),
528 m_vtx0Daug1MassHypo(-1),
529 m_vtx1Daug1MassHypo(-1),
530 m_vtx1Daug2MassHypo(-1),
531 m_vtx1Daug3MassHypo(-1),
535 m_iVertexFitter(
"Trk::TrkVKalVrtFitter"),
536 m_pvRefitter(
"Analysis::PrimaryVertexRefitter",
this),
537 m_V0Tools(
"Trk::V0Tools"),
538 m_CascadeTools(
"DerivationFramework::CascadeTools"),
539 m_muonCollectionKey(
"StacoMuonCollection"),
540 m_trkSelector(
"InDet::TrackSelectorTool"),
544 m_useCombMeasurement(
false)
546 declareProperty(
"DxVertices", m_vertexDxContainerKey);
547 declareProperty(
"VxPrimaryCandidateName", m_VxPrimaryCandidateName);
548 declareProperty(
"RefPVContainerName", m_refPVContainerName =
"RefittedPrimaryVertices");
549 declareProperty(
"DxMassLowerCut", m_DxMassLower);
550 declareProperty(
"DxMassUpperCut", m_DxMassUpper);
551 declareProperty(
"MassLowerCut", m_MassLower);
552 declareProperty(
"MassUpperCut", m_MassUpper);
553 declareProperty(
"HypothesisName", m_hypoName =
"B");
554 declareProperty(
"Vtx0MassHypo", m_vtx0MassHypo);
555 declareProperty(
"Vtx1MassHypo", m_vtx1MassHypo);
556 declareProperty(
"Vtx0Daug1MassHypo", m_vtx0Daug1MassHypo);
557 declareProperty(
"Vtx1Daug1MassHypo", m_vtx1Daug1MassHypo);
558 declareProperty(
"Vtx1Daug2MassHypo", m_vtx1Daug2MassHypo);
559 declareProperty(
"Vtx1Daug3MassHypo", m_vtx1Daug3MassHypo);
560 declareProperty(
"DxHypothesis", m_Dx_pid);
561 declareProperty(
"ApplyDxMassConstraint", m_constrDx);
562 declareProperty(
"Chi2Cut", m_chi2cut);
563 declareProperty(
"RefitPV", m_refitPV =
true);
564 declareProperty(
"MaxnPV", m_PV_max = 999);
565 declareProperty(
"MinNTracksInPV", m_PV_minNTracks = 0);
566 declareProperty(
"DoVertexType", m_DoVertexType = 7);
567 declareProperty(
"TrkVertexFitterTool", m_iVertexFitter);
568 declareProperty(
"PVRefitter", m_pvRefitter);
569 declareProperty(
"V0Tools", m_V0Tools);
570 declareProperty(
"CascadeTools", m_CascadeTools);
571 declareProperty(
"CascadeVertexCollections", m_cascadeOutputsKeys);
572 declareProperty(
"muonCollectionKey", m_muonCollectionKey);
573 declareProperty(
"TrackSelectorTool", m_trkSelector);
574 declareProperty(
"muonThresholdPt", m_thresholdPt);
575 declareProperty(
"useMCPCuts", m_mcpCuts);
576 declareProperty(
"combOnly", m_combOnly);
577 declareProperty(
"useCombinedMeasurement", m_useCombMeasurement);
585 assert(cascadeinfoContainer!=
nullptr);
596 std::vector<const xAOD::TrackParticle*> tracksMu;
597 std::vector<const xAOD::TrackParticle*> tracksDx;
598 std::vector<double> massesMu;
600 std::vector<double> massesDx;
604 std::vector<double> massesDm;
608 std::vector<double> Masses;
620 return StatusCode::SUCCESS;;
627 typedef std::vector<const xAOD::Muon*>
MuonBag;
631 MuonBag theMuonsAfterSelection;
632 for (
auto mu : *importedMuonCollection) {
633 if (
mu ==
nullptr )
continue;
634 if (!
mu->inDetTrackParticleLink().isValid())
continue;
636 if ( muonTrk==
nullptr)
continue;
643 theMuonsAfterSelection.push_back(
mu);
645 if (theMuonsAfterSelection.size() == 0)
return StatusCode::SUCCESS;;
646 ATH_MSG_DEBUG(
"Number of muons after selection: " << theMuonsAfterSelection.size());
651 std::vector<const xAOD::Vertex*> selectedDxCandidates;
654 for(
auto vxcItr : *dxContainer){
661 if(!flagAcc1(*vtx))
continue;
671 if(!flagAcc1(*vtx)) isDp =
false;
674 if(!flagAcc2(*vtx)) isDm =
false;
676 if(!(isDp||isDm))
continue;
681 ATH_MSG_DEBUG(
" Original Dx/Lambda_c candidate rejected by the track's cut level - loose ");
685 ATH_MSG_DEBUG(
" Original Dx/Lambda_c candidate rejected by the track's cut level - loose ");
689 ATH_MSG_DEBUG(
" Original Dx/Lambda_c candidate rejected by the track's cut level - loose ");
694 if( std::abs( vxcItr->trackParticle(0)->charge()+vxcItr->trackParticle(1)->charge()+vxcItr->trackParticle(2)->charge() ) != 1 ){
695 ATH_MSG_DEBUG(
" Original Dx/Lambda_c candidate rejected by the charge requirement: "
696 << vxcItr->trackParticle(0)->charge() <<
", " << vxcItr->trackParticle(1)->charge() <<
", " << vxcItr->trackParticle(2)->charge() );
702 if( (std::abs(
m_Dx_pid) == 411 || std::abs(
m_Dx_pid) == 4122) && vxcItr->trackParticle(2)->charge()<0)
703 mass_D =
m_V0Tools->invariantMass(vxcItr,massesDm);
705 mass_D =
m_V0Tools->invariantMass(vxcItr,massesDx);
708 ATH_MSG_DEBUG(
" Original D_(s)/Lambda_c candidate rejected by the mass cut: mass = "
715 TLorentzVector p4Kp_in, p4Km_in;
716 p4Kp_in.SetPtEtaPhiM( vxcItr->trackParticle(0)->pt(),
717 vxcItr->trackParticle(0)->eta(),
719 p4Km_in.SetPtEtaPhiM( vxcItr->trackParticle(1)->pt(),
720 vxcItr->trackParticle(1)->eta(),
722 double mass_phi = (p4Kp_in + p4Km_in).M();
724 if(mass_phi > 1600) {
725 ATH_MSG_DEBUG(
" Original phi candidate rejected by the mass cut: mass = " << mass_phi );
733 if (vxcItr->trackParticle(2)->pt() < 2400)
continue;
737 if(std::abs(
m_Dx_pid)==411 && selectedDxCandidates.size()>0) {
739 for(
auto dxItr : selectedDxCandidates){
740 if(vxcItr->trackParticle(2)->charge()>0) {
741 if ( vxcItr->trackParticle(2) == dxItr->trackParticle(0) && vxcItr->trackParticle(0) == dxItr->trackParticle(2)
742 && vxcItr->trackParticle(1) == dxItr->trackParticle(1) ) Dpmcopy =
true;
744 if ( vxcItr->trackParticle(2) == dxItr->trackParticle(1) && vxcItr->trackParticle(1) == dxItr->trackParticle(2)
745 && vxcItr->trackParticle(0) == dxItr->trackParticle(0) ) Dpmcopy =
true;
748 if (Dpmcopy)
continue;
751 selectedDxCandidates.push_back(vxcItr);
753 if(selectedDxCandidates.size()<1)
return StatusCode::SUCCESS;
759 for(
auto muItr : theMuonsAfterSelection){
762 auto TrkMuon = muItr->trackParticle( xAOD::Muon::InnerDetectorTrackParticle );
763 tracksMu.push_back(TrkMuon);
764 if (tracksMu.size() != 1 || massesMu.size() != 1 ) {
769 for(
auto dxItr : selectedDxCandidates){
771 if(
std::find(tracksMu.cbegin(), tracksMu.cend(), dxItr->trackParticle(0)) != tracksMu.cend())
continue;
772 if(
std::find(tracksMu.cbegin(), tracksMu.cend(), dxItr->trackParticle(1)) != tracksMu.cend())
continue;
773 if(
std::find(tracksMu.cbegin(), tracksMu.cend(), dxItr->trackParticle(2)) != tracksMu.cend())
continue;
775 size_t dxTrkNum = dxItr->nTrackParticles();
777 for(
unsigned int it=0;
it<dxTrkNum;
it++) tracksDx.push_back(dxItr->trackParticle(
it));
778 if (tracksDx.size() != 3 || massesDx.size() != 3 ) {
782 if(
std::find(tracksDx.cbegin(), tracksDx.cend(), TrkMuon) != tracksDx.cend())
continue;
784 ATH_MSG_DEBUG(
"Using tracks" << tracksMu[0] <<
", " << tracksDx[0] <<
", " << tracksDx[1] <<
", " << tracksDx[2]);
793 std::vector<Trk::VertexID> vrtList;
797 if((std::abs(
m_Dx_pid)==411 || std::abs(
m_Dx_pid)==4122 ) && dxItr->trackParticle(2)->charge()<0)
802 if((std::abs(
m_Dx_pid)==411 || std::abs(
m_Dx_pid)==4122 ) && dxItr->trackParticle(2)->charge()<0)
807 vrtList.push_back(vID);
821 << ((
result->vertices())[0])->trackParticle(1) <<
", "
822 << ((
result->vertices())[0])->trackParticle(2) <<
", "
823 << ((
result->vertices())[1])->trackParticle(0));
826 result->setSVOwnership(
true);
833 const std::vector<xAOD::Vertex*> &cascadeVertices =
result->vertices();
834 const std::vector< std::vector<TLorentzVector> > &moms =
result->getParticleMoms();
844 if (pvContainer->
size()==0){
846 return StatusCode::RECOVERABLE;
848 primaryVertex = (*pvContainer)[0];
853 TLorentzVector p4Kp_in, p4Km_in;
854 p4Kp_in.SetPtEtaPhiM( cascadeVertices[0]->trackParticle(0)->
pt(),
855 cascadeVertices[0]->trackParticle(0)->eta(),
857 p4Km_in.SetPtEtaPhiM( cascadeVertices[0]->trackParticle(1)->
pt(),
858 cascadeVertices[0]->trackParticle(1)->eta(),
860 double mass_phi = (p4Kp_in + p4Km_in).M();
862 if(mass_phi > 1100) {
863 ATH_MSG_DEBUG(
" Original phi candidate rejected by the mass cut: mass = " << mass_phi );
871 if (
m_CascadeTools->lxy(moms[1],cascadeVertices[1],primaryVertex) > 0.09){
873 if (
m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]) > 0){
874 if (std::abs(
m_Dx_pid)!=411) cascadeinfoContainer->push_back(
result.release());
875 else if (
m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]) > 0.09){
876 cascadeinfoContainer->push_back(
result.release());
893 ATH_MSG_DEBUG(
"cascadeinfoContainer size " << cascadeinfoContainer->size());
895 return StatusCode::SUCCESS;