21 #include "AthLinks/ElementLink.h"
30 constexpr
double muMass = 105.658;
31 constexpr
double kMass = 493.677;
32 constexpr
double piMass = 139.57;
33 constexpr
double pMass = 938.272;
52 ATH_MSG_FATAL(
"Invalid number of elements given for manualMass hypothesis - needs 4");
53 return StatusCode::FAILURE;
61 return StatusCode::FAILURE;
65 return StatusCode::FAILURE;
83 return StatusCode::SUCCESS;
92 m_oppChargesOnly(true),
93 m_sameChargesOnly(false),
94 m_trkThresholdPt(0.0),
99 m_jpsiCollectionKey(
"JpsiCandidates"),
100 m_jpsiMassUpper(0.0),
101 m_jpsiMassLower(0.0),
102 m_TrkParticleCollection(
"InDetTrackParticles"),
103 m_TrkParticleGSFCollection(
"GSFTrackParticles"),
105 m_excludeJpsiMuonsOnly(true),
106 m_excludeCrossJpsiTracks(false),
107 m_iVertexFitter(
"Trk::TrkVKalVrtFitter"),
108 m_trkSelector(
"InDet::TrackSelectorTool"),
109 m_useMassConst(true),
110 m_altMassConst(-1.0),
111 m_diTrackMassUpper(-1.0),
112 m_diTrackMassLower(-1.0),
115 m_trkQuadrupletMassUpper(-1.0),
116 m_trkQuadrupletMassLower(-1.0),
117 m_trkQuadrupletPt(-1.0),
118 m_finalDiTrackMassUpper(-1.0),
119 m_finalDiTrackMassLower(-1.0),
120 m_finalDiTrackPt(-1.0),
123 m_candidateLimit(std::numeric_limits<size_t>::
max())
125 declareInterface<JpsiPlus2Tracks>(
this);
183 return StatusCode::FAILURE;
185 importedJpsiCollection = jpsihandle.
cptr();
188 ATH_MSG_DEBUG(
"VxCandidate container size " << importedJpsiCollection->size());
195 return StatusCode::FAILURE;
197 importedTrackCollection = trackshandle.
cptr();
200 ATH_MSG_DEBUG(
"Track container size "<< importedTrackCollection->size());
206 importedGSFTrackCollection =
h.cptr();
215 return StatusCode::FAILURE;
217 importedMuonCollection =
h.cptr();
220 ATH_MSG_DEBUG(
"Muon container size "<< importedMuonCollection->size());
224 typedef std::vector<const xAOD::TrackParticle*>
TrackBag;
228 for (
auto trkPBItr=importedTrackCollection->cbegin(); trkPBItr!=importedTrackCollection->cend(); ++trkPBItr) {
235 if (
m_trkSelector->decision(*
tp, NULL) ) theIDTracksAfterSelection.push_back(
tp);
237 if (theIDTracksAfterSelection.size() == 0)
return StatusCode::SUCCESS;
238 ATH_MSG_DEBUG(
"Number of tracks after ID trkSelector: " << theIDTracksAfterSelection.size());
241 std::vector<const xAOD::Vertex*> selectedJpsiCandidates;
242 std::vector<const xAOD::TrackParticle*> jpsiTracks;
243 for(
auto vxcItr=importedJpsiCollection->cbegin(); vxcItr!=importedJpsiCollection->cend(); ++vxcItr) {
252 selectedJpsiCandidates.push_back(*vxcItr);
259 jpsiTracks.push_back(jpsiTP1);
260 jpsiTracks.push_back(jpsiTP2);
263 ATH_MSG_DEBUG(
"selectedJpsiCandidates: " << selectedJpsiCandidates.size());
270 std::vector<const xAOD::TrackParticle*> QuadletTracks(4,
nullptr);
272 std::vector<double> massCuts;
276 for(
auto muon : *importedMuonCollection){
277 if(!
muon->inDetTrackParticleLink().isValid())
continue;
278 auto track =
muon->trackParticle( xAOD::Muon::InnerDetectorTrackParticle );
279 if(
track==
nullptr)
continue;
281 muonTracks.push_back(
track);
285 for(
auto jpsiItr=selectedJpsiCandidates.begin(); jpsiItr!=selectedJpsiCandidates.end(); ++jpsiItr) {
290 QuadletTracks[0] = jpsiTP1;
291 QuadletTracks[1] = jpsiTP2;
295 jpsiTracks.resize(2);
296 jpsiTracks[0] = jpsiTP1;
297 jpsiTracks[1] = jpsiTP2;
301 for (
TrackBag::iterator trkItr1=theIDTracksAfterSelection.begin(); trkItr1<theIDTracksAfterSelection.end(); ++trkItr1) {
303 int linkedMuonTrk1 = 0;
306 if (linkedMuonTrk1)
ATH_MSG_DEBUG(
"This id track 1 is muon track!");
309 if (linkedMuonTrk1)
ATH_MSG_DEBUG(
"ID track 1 removed: id track is selected to build Jpsi!");
316 std::abs((*trkItr1)->z0() + (*trkItr1)->vz() - (*jpsiItr)->z()) >
m_trkDeltaZ )
319 for (
TrackBag::iterator trkItr2=trkItr1+1; trkItr2!=theIDTracksAfterSelection.end(); ++trkItr2) {
322 return StatusCode::SUCCESS;
328 if (linkedMuonTrk2)
ATH_MSG_DEBUG(
"This id track 2 is muon track!");
330 if (linkedMuonTrk2)
ATH_MSG_DEBUG(
"ID track 2 removed: id track is selected to build Jpsi Vtx!");
341 std::abs((*trkItr2)->z0() + (*trkItr2)->vz() - (*jpsiItr)->z()) >
m_trkDeltaZ )
350 if((*trkItr2)->qOverP()>0) {
351 QuadletTracks[2] = *trkItr2;
352 QuadletTracks[3] = *trkItr1;
354 QuadletTracks[2] = *trkItr1;
355 QuadletTracks[3] = *trkItr2;
361 bool passesDiTrack(
true);
378 if (!passesDiTrack)
continue;
381 bool passes4TrackMass(
true);
399 if (!passes4TrackMass)
continue;
402 std::unique_ptr<xAOD::Vertex> bVertex (
fit(QuadletTracks, importedTrackCollection, importedGSFTrackCollection));
403 if(!bVertex)
continue;
407 if(!chi2CutPassed) {
ATH_MSG_DEBUG(
"Chi Cut failed!");
continue; }
411 if(!passesCuts)
continue;
414 std::vector<const xAOD::Vertex*> theJpsiPreceding;
415 theJpsiPreceding.push_back(*jpsiItr);
423 return StatusCode::SUCCESS;
449 constexpr
double jpsiTableMass = 3096.916;
451 std::array<int,2>
indices= {1, 2};
467 std::vector<ElementLink<DataVector<xAOD::TrackParticle> > > newLinkVector;
472 newLinkVector.push_back( mylink );
488 double qOverP1=trk1->
qOverP();
489 double qOverP2=trk2->
qOverP();
490 bool opposite = (qOverP1*qOverP2<0.0);
504 const auto trk1V = trk1->
p4();
505 double px1 = trk1V.Px();
506 double py1 = trk1V.Py();
507 double pz1 = trk1V.Pz();
508 double e1 = sqrt(px1*px1+py1*py1+pz1*pz1+mass1*mass1);
509 const auto trk2V = trk2->
p4();
510 double px2 = trk2V.Px();
511 double py2 = trk2V.Py();
512 double pz2 = trk2V.Pz();
513 double e2 = sqrt(px2*px2+py2*py2+pz2*pz2+mass2*mass2);
514 double pxSum=px1+px2;
515 double pySum=py1+py2;
516 double pzSum=pz1+pz2;
518 double M=sqrt((eSum*eSum)-(pxSum*pxSum)-(pySum*pySum)-(pzSum*pzSum));
525 const std::vector<double> &
masses)
527 assert(trk.size() ==
masses.size() && trk.size()==4);
528 const auto trk1V = trk[0]->p4();
529 double px1 = trk1V.Px();
530 double py1 = trk1V.Py();
531 double pz1 = trk1V.Pz();
534 const auto trk2V = trk[1]->p4();
535 double px2 = trk2V.Px();
536 double py2 = trk2V.Py();
537 double pz2 = trk2V.Pz();
540 const auto trk3V = trk[2]->p4();
541 double px3 = trk3V.Px();
542 double py3 = trk3V.Py();
543 double pz3 = trk3V.Pz();
546 const auto trk4V = trk[3]->p4();
547 double px4 = trk4V.Px();
548 double py4 = trk4V.Py();
549 double pz4 = trk4V.Pz();
552 double pxSum=px1+px2+px3+px4;
553 double pySum=py1+py2+py3+py4;
554 double pzSum=pz1+pz2+pz3+pz4;
557 double M=sqrt((eSum*eSum)-(pxSum*pxSum)-(pySum*pySum)-(pzSum*pzSum));
567 double bPt = bMomentum.Pt();
569 double bMass = bMomentum.M();
570 ATH_MSG_DEBUG(
"Candidate pt/mass under " <<
str <<
" track mass hypothesis is " << bPt <<
" / " << bMass);
576 double bDiTrkPt = (tr1+tr2).
Pt();
577 double bDiTrkMass = (tr1+tr2).M();