21 #include "AthLinks/ElementLink.h"
53 ATH_MSG_FATAL(
"Invalid number of elements given for manualMass hypothesis - needs 4");
54 return StatusCode::FAILURE;
62 return StatusCode::FAILURE;
66 return StatusCode::FAILURE;
84 return StatusCode::SUCCESS;
93 m_oppChargesOnly(true),
94 m_sameChargesOnly(false),
95 m_trkThresholdPt(0.0),
100 m_jpsiCollectionKey(
"JpsiCandidates"),
101 m_jpsiMassUpper(0.0),
102 m_jpsiMassLower(0.0),
103 m_TrkParticleCollection(
"InDetTrackParticles"),
104 m_TrkParticleGSFCollection(
"GSFTrackParticles"),
106 m_excludeJpsiMuonsOnly(true),
107 m_excludeCrossJpsiTracks(false),
108 m_iVertexFitter(
"Trk::TrkVKalVrtFitter"),
109 m_trkSelector(
"InDet::TrackSelectorTool"),
110 m_useMassConst(true),
111 m_altMassConst(-1.0),
112 m_diTrackMassUpper(-1.0),
113 m_diTrackMassLower(-1.0),
116 m_trkQuadrupletMassUpper(-1.0),
117 m_trkQuadrupletMassLower(-1.0),
118 m_trkQuadrupletPt(-1.0),
119 m_finalDiTrackMassUpper(-1.0),
120 m_finalDiTrackMassLower(-1.0),
121 m_finalDiTrackPt(-1.0),
124 m_candidateLimit(std::numeric_limits<size_t>::
max())
126 declareInterface<JpsiPlus2Tracks>(
this);
184 return StatusCode::FAILURE;
186 importedJpsiCollection = jpsihandle.
cptr();
189 ATH_MSG_DEBUG(
"VxCandidate container size " << importedJpsiCollection->size());
196 return StatusCode::FAILURE;
198 importedTrackCollection = trackshandle.
cptr();
201 ATH_MSG_DEBUG(
"Track container size "<< importedTrackCollection->size());
207 importedGSFTrackCollection =
h.cptr();
216 return StatusCode::FAILURE;
218 importedMuonCollection =
h.cptr();
221 ATH_MSG_DEBUG(
"Muon container size "<< importedMuonCollection->size());
225 typedef std::vector<const xAOD::TrackParticle*>
TrackBag;
229 for (
auto trkPBItr=importedTrackCollection->cbegin(); trkPBItr!=importedTrackCollection->cend(); ++trkPBItr) {
236 if (
m_trkSelector->decision(*
tp, NULL) ) theIDTracksAfterSelection.push_back(
tp);
238 if (theIDTracksAfterSelection.size() == 0)
return StatusCode::SUCCESS;
239 ATH_MSG_DEBUG(
"Number of tracks after ID trkSelector: " << theIDTracksAfterSelection.size());
242 std::vector<const xAOD::Vertex*> selectedJpsiCandidates;
243 std::vector<const xAOD::TrackParticle*> jpsiTracks;
244 for(
auto vxcItr=importedJpsiCollection->cbegin(); vxcItr!=importedJpsiCollection->cend(); ++vxcItr) {
253 selectedJpsiCandidates.push_back(*vxcItr);
260 jpsiTracks.push_back(jpsiTP1);
261 jpsiTracks.push_back(jpsiTP2);
264 ATH_MSG_DEBUG(
"selectedJpsiCandidates: " << selectedJpsiCandidates.size());
271 std::vector<const xAOD::TrackParticle*> QuadletTracks(4,
nullptr);
273 std::vector<double> massCuts;
277 for(
auto muon : *importedMuonCollection){
278 if(!
muon->inDetTrackParticleLink().isValid())
continue;
279 auto track =
muon->trackParticle( xAOD::Muon::InnerDetectorTrackParticle );
280 if(
track==
nullptr)
continue;
282 muonTracks.push_back(
track);
286 for(
auto jpsiItr=selectedJpsiCandidates.begin(); jpsiItr!=selectedJpsiCandidates.end(); ++jpsiItr) {
291 QuadletTracks[0] = jpsiTP1;
292 QuadletTracks[1] = jpsiTP2;
296 jpsiTracks.resize(2);
297 jpsiTracks[0] = jpsiTP1;
298 jpsiTracks[1] = jpsiTP2;
302 for (
TrackBag::iterator trkItr1=theIDTracksAfterSelection.begin(); trkItr1<theIDTracksAfterSelection.end(); ++trkItr1) {
304 int linkedMuonTrk1 = 0;
307 if (linkedMuonTrk1)
ATH_MSG_DEBUG(
"This id track 1 is muon track!");
310 if (linkedMuonTrk1)
ATH_MSG_DEBUG(
"ID track 1 removed: id track is selected to build Jpsi!");
317 std::abs((*trkItr1)->z0() + (*trkItr1)->vz() - (*jpsiItr)->z()) >
m_trkDeltaZ )
320 for (
TrackBag::iterator trkItr2=trkItr1+1; trkItr2!=theIDTracksAfterSelection.end(); ++trkItr2) {
323 return StatusCode::SUCCESS;
329 if (linkedMuonTrk2)
ATH_MSG_DEBUG(
"This id track 2 is muon track!");
331 if (linkedMuonTrk2)
ATH_MSG_DEBUG(
"ID track 2 removed: id track is selected to build Jpsi Vtx!");
342 std::abs((*trkItr2)->z0() + (*trkItr2)->vz() - (*jpsiItr)->z()) >
m_trkDeltaZ )
351 if((*trkItr2)->qOverP()>0) {
352 QuadletTracks[2] = *trkItr2;
353 QuadletTracks[3] = *trkItr1;
355 QuadletTracks[2] = *trkItr1;
356 QuadletTracks[3] = *trkItr2;
362 bool passesDiTrack(
true);
379 if (!passesDiTrack)
continue;
382 bool passes4TrackMass(
true);
400 if (!passes4TrackMass)
continue;
403 std::unique_ptr<xAOD::Vertex> bVertex (
fit(QuadletTracks, importedTrackCollection, importedGSFTrackCollection));
404 if(!bVertex)
continue;
408 if(!chi2CutPassed) {
ATH_MSG_DEBUG(
"Chi Cut failed!");
continue; }
412 if(!passesCuts)
continue;
415 std::vector<const xAOD::Vertex*> theJpsiPreceding;
416 theJpsiPreceding.push_back(*jpsiItr);
424 return StatusCode::SUCCESS;
452 std::array<int,2>
indices= {1, 2};
468 std::vector<ElementLink<DataVector<xAOD::TrackParticle> > > newLinkVector;
473 newLinkVector.push_back( mylink );
489 double qOverP1=trk1->
qOverP();
490 double qOverP2=trk2->
qOverP();
491 bool opposite = (qOverP1*qOverP2<0.0);
505 const auto trk1V = trk1->
p4();
506 double px1 = trk1V.Px();
507 double py1 = trk1V.Py();
508 double pz1 = trk1V.Pz();
509 double e1 = sqrt(px1*px1+py1*py1+pz1*pz1+mass1*mass1);
510 const auto trk2V = trk2->
p4();
511 double px2 = trk2V.Px();
512 double py2 = trk2V.Py();
513 double pz2 = trk2V.Pz();
514 double e2 = sqrt(px2*px2+py2*py2+pz2*pz2+mass2*mass2);
515 double pxSum=px1+px2;
516 double pySum=py1+py2;
517 double pzSum=pz1+pz2;
519 double M=sqrt((eSum*eSum)-(pxSum*pxSum)-(pySum*pySum)-(pzSum*pzSum));
526 const std::vector<double> &
masses)
528 assert(trk.size() ==
masses.size() && trk.size()==4);
529 const auto trk1V = trk[0]->p4();
530 double px1 = trk1V.Px();
531 double py1 = trk1V.Py();
532 double pz1 = trk1V.Pz();
535 const auto trk2V = trk[1]->p4();
536 double px2 = trk2V.Px();
537 double py2 = trk2V.Py();
538 double pz2 = trk2V.Pz();
541 const auto trk3V = trk[2]->p4();
542 double px3 = trk3V.Px();
543 double py3 = trk3V.Py();
544 double pz3 = trk3V.Pz();
547 const auto trk4V = trk[3]->p4();
548 double px4 = trk4V.Px();
549 double py4 = trk4V.Py();
550 double pz4 = trk4V.Pz();
553 double pxSum=px1+px2+px3+px4;
554 double pySum=py1+py2+py3+py4;
555 double pzSum=pz1+pz2+pz3+pz4;
558 double M=sqrt((eSum*eSum)-(pxSum*pxSum)-(pySum*pySum)-(pzSum*pzSum));
568 double bPt = bMomentum.Pt();
570 double bMass = bMomentum.M();
571 ATH_MSG_DEBUG(
"Candidate pt/mass under " <<
str <<
" track mass hypothesis is " << bPt <<
" / " << bMass);
577 double bDiTrkPt = (tr1+tr2).
Pt();
578 double bDiTrkMass = (tr1+tr2).M();