22 #include "GaudiKernel/IPartPropSvc.h"
23 #include "AthLinks/ElementLink.h"
24 #include "HepPDT/ParticleDataTable.hh"
52 IPartPropSvc* partPropSvc = 0;
53 StatusCode sc = service(
"PartPropSvc", partPropSvc,
true);
55 ATH_MSG_ERROR(
"Could not initialize Particle Properties Service");
56 return StatusCode::FAILURE;
58 auto particleDataTable = partPropSvc->PDT();
59 const HepPDT::ParticleData* pd_el = particleDataTable->particle(MC::ELECTRON);
63 if (
m_doTagAndProbe)
ATH_MSG_WARNING(
"You have requested tag and probe mode. Duplicate mu+trk pairs WILL be allowed, charge ordering WILL NOT be done. Tag track will be first in each candidate");
67 bool illogicalOptions(
false);
69 ATH_MSG_WARNING(
"You are requesting incompatible combinations of muons and tracks in the pairs. JpsiEECandidates will be EMPTY!");
70 illogicalOptions=
true;
73 ATH_MSG_WARNING(
"You are requesting Tag and Probe analysis but have not requested mu+trk mode. This is impossible. JpsiEECandidates will be EMPTY!");
74 illogicalOptions=
true;
77 ATH_MSG_WARNING(
"You are requesting incompatible combinations of charges in the pairs. JpsiEECandidates will be EMPTY!");
78 illogicalOptions=
true;
81 ATH_MSG_WARNING(
"You are requesting same-sign or all-sign combinations in a tag and probe analysis. This doesn't make sense. JpsiEECandidates will be EMPTY!");
82 illogicalOptions=
true;
84 if (illogicalOptions)
return StatusCode::FAILURE;
90 return StatusCode::SUCCESS;
99 m_allElectrons(false),
100 m_useTrackMeasurement(true),
101 m_useV0Fitter(false),
107 m_trkThresholdPt(0.0),
108 m_invMassUpper(100000.0),
110 m_collAngleTheta(0.0),
116 m_electronCollectionKey(
"Electrons"),
117 m_TrkParticleCollection(
"InDetTrackParticles"),
118 m_iVertexFitter(
"Trk::TrkVKalVrtFitter"),
119 m_iV0VertexFitter(
"Trk::V0VertexFitter"),
120 m_trkSelector(
"InDet::TrackSelectorTool"),
121 m_vertexEstimator(
"InDet::VertexPointEstimator"),
123 m_elSelection(
"d0_or_nod0"),
124 m_doTagAndProbe(false)
127 declareInterface<JpsiFinder_ee>(
this);
177 return StatusCode::SUCCESS;;
179 importedElectronCollection = ehandle.
cptr();
189 return StatusCode::SUCCESS;;
191 importedTrackCollection = thandle.
cptr();
195 typedef std::vector<const xAOD::TrackParticle*>
TrackBag;
196 typedef std::vector<const xAOD::Electron*> ElectronBag;
203 for (trkCItr=importedTrackCollection->
begin(); trkCItr!=importedTrackCollection->
end(); ++trkCItr) {
207 theIDTracksAfterSelection.push_back(TP);
209 if (theIDTracksAfterSelection.size() == 0)
return StatusCode::SUCCESS;;
210 ATH_MSG_DEBUG(
"Number of tracks after ID track selection: " << theIDTracksAfterSelection.size());
214 ElectronBag theElectronsAfterSelection;
217 for (elItr=importedElectronCollection->
begin(); elItr!=importedElectronCollection->
end(); ++elItr) {
218 if ( *elItr == NULL )
continue;
219 if (!(*elItr)->trackParticleLink().isValid())
continue;
221 elTrk = (*elItr)->trackParticleLink().cachedElement();
223 if ( elTrk==NULL)
continue;
228 theElectronsAfterSelection.push_back(*elItr);
230 if (theElectronsAfterSelection.size() == 0)
return StatusCode::SUCCESS;;
231 ATH_MSG_DEBUG(
"Number of electrons after selection: " << theElectronsAfterSelection.size());
235 std::vector<JpsiEECandidate> jpsiCandidates;
236 if (
m_elel) jpsiCandidates =
getPairs(theElectronsAfterSelection);
240 ATH_MSG_DEBUG(
"Number of pairs with ee from a B decay: " << jpsiCandidates.size() );
252 for (jpsiItr=jpsiCandidates.begin(); jpsiItr!=jpsiCandidates.end(); ++jpsiItr) {
253 (*jpsiItr).collection1 = importedTrackCollection;
254 (*jpsiItr).collection2 = importedTrackCollection;
259 for (jpsiItr=jpsiCandidates.begin(); jpsiItr!=jpsiCandidates.end(); ++jpsiItr) {
261 (*jpsiItr).trackParticle1 = (*jpsiItr).el1->trackParticleLink().cachedElement();
262 (*jpsiItr).trackParticle2 = (*jpsiItr).el2->trackParticleLink().cachedElement();
263 (*jpsiItr).collection1 = importedTrackCollection;
264 (*jpsiItr).collection2 = importedTrackCollection;
266 ATH_MSG_WARNING(
"Not setup for non-track electron measurements yet....");
275 std::vector<int> listToDelete;
276 std::vector<int>::reverse_iterator ii;
277 for(jpsiItr=jpsiCandidates.begin(); jpsiItr!=jpsiCandidates.end();++jpsiItr,++
index) {
278 if( (fabs((*jpsiItr).trackParticle1->pt()) <
m_higherPt) && (fabs((*jpsiItr).trackParticle2->pt()) <
m_higherPt) ) listToDelete.push_back(
index);
280 for (ii=listToDelete.rbegin(); ii!=listToDelete.rend(); ++ii) {
281 jpsiCandidates.erase(jpsiCandidates.begin() + (*ii) );
283 ATH_MSG_DEBUG(
"Number of candidates after higherPt cut: " << jpsiCandidates.size() );
287 std::vector<JpsiEECandidate> sortedJpsiEECandidates;
291 ATH_MSG_DEBUG(
"Number of candidates after charge selection: " << sortedJpsiEECandidates.size() );
296 std::vector<int> listToDelete;
297 std::vector<int>::reverse_iterator ii;
298 for(jpsiItr=sortedJpsiEECandidates.begin(); jpsiItr!=sortedJpsiEECandidates.end();++jpsiItr,++
index) {
299 double deltatheta = fabs( (*jpsiItr).trackParticle1->theta() - (*jpsiItr).trackParticle2->theta() );
301 double deltaphi = std::abs(
xAOD::P4Helpers::deltaPhi((*jpsiItr).trackParticle1->phi0() , (*jpsiItr).trackParticle2->phi0()));
305 for (ii=listToDelete.rbegin(); ii!=listToDelete.rend(); ++ii) {
306 sortedJpsiEECandidates.erase(sortedJpsiEECandidates.begin() + (*ii) );
308 ATH_MSG_DEBUG(
"Number of collimated candidates: " << sortedJpsiEECandidates.size() );
312 std::vector<double> trkMasses;
317 std::vector<int> listToDelete;
318 std::vector<int>::reverse_iterator ii;
319 for(jpsiItr=sortedJpsiEECandidates.begin(); jpsiItr!=sortedJpsiEECandidates.end(); ++jpsiItr,++
index) {
323 listToDelete.push_back(
index);
326 for (ii=listToDelete.rbegin(); ii!=listToDelete.rend(); ++ii) {
327 sortedJpsiEECandidates.erase(sortedJpsiEECandidates.begin() + (*ii) );
329 ATH_MSG_DEBUG(
"Number of candidates passing invariant mass selection: " << sortedJpsiEECandidates.size() );
332 if (sortedJpsiEECandidates.size() == 0)
return StatusCode::SUCCESS;;
336 for(jpsiItr=sortedJpsiEECandidates.begin(); jpsiItr!=sortedJpsiEECandidates.end(); ++jpsiItr) {
338 std::vector<const xAOD::TrackParticle*> theTracks; theTracks.clear();
339 theTracks.push_back((*jpsiItr).trackParticle1);
340 theTracks.push_back((*jpsiItr).trackParticle2);
341 ATH_MSG_DEBUG(
"theTracks size (should be two!) " << theTracks.size() <<
" being vertexed with tracks " << importedTrackCollection);
342 std::unique_ptr<xAOD::Vertex> myVxCandidate{
fit(theTracks,importedTrackCollection)};
343 if (myVxCandidate != 0) {
345 double chi2 = myVxCandidate->chiSquared();
352 std::vector<const xAOD::Electron*> theStoredElectrons;
353 theStoredElectrons.push_back((*jpsiItr).el1);
354 if (
m_elel) theStoredElectrons.push_back((*jpsiItr).el2);
355 jpsiHelper.
setElectrons(theStoredElectrons,importedElectronCollection);
358 vxContainer.
push_back(std::move(myVxCandidate));
369 return StatusCode::SUCCESS;;
385 if(concreteVertexFitter == 0) {
386 ATH_MSG_FATAL(
"The vertex fitter passed is not a V0 Vertex Fitter");
396 if (errorcode != 0) {startingPoint(0) = 0.0; startingPoint(1) = 0.0; startingPoint(2) = 0.0;}
401 if(myVxCandidate != 0){
402 std::vector<ElementLink<DataVector<xAOD::TrackParticle> > > newLinkVector;
406 newLinkVector.push_back( mylink ); }
414 return myVxCandidate;
417 ATH_MSG_DEBUG(
"Initial fit was a success! " << myVxCandidate);
419 if(myVxCandidate != 0){
420 std::vector<ElementLink<DataVector<xAOD::TrackParticle> > > newLinkVector;
424 newLinkVector.push_back( mylink );
432 return myVxCandidate;
450 std::vector<JpsiEECandidate> myPairs;
452 std::vector<const xAOD::TrackParticle*>::const_iterator outerItr;
453 std::vector<const xAOD::TrackParticle*>::const_iterator innerItr;
455 if(TracksIn.size()>=2){
456 for(outerItr=TracksIn.begin();outerItr<TracksIn.end();++outerItr){
457 for(innerItr=(outerItr+1);innerItr!=TracksIn.end();++innerItr){
461 myPairs.push_back(pair);
477 std::vector<JpsiEECandidate> myPairs;
479 std::vector<const xAOD::Electron*>::const_iterator outerItr;
480 std::vector<const xAOD::Electron*>::const_iterator innerItr;
482 if(electronsIn.size()>=2){
483 for(outerItr=electronsIn.begin();outerItr<electronsIn.end();++outerItr){
484 for(innerItr=(outerItr+1);innerItr!=electronsIn.end();++innerItr){
487 pair.
el1 = *innerItr;
488 pair.
el2 = *outerItr;
490 myPairs.push_back(pair);
506 std::vector<JpsiEECandidate> myPairs;
510 std::vector<const xAOD::TrackParticle*> tracksToKeep;
512 if(tracks.size()>=1 &&
electrons.size()>=1){
514 bool trackIsElectron(
false);
516 if ( ele->trackParticleLink().cachedElement() == trk ) {
517 trackIsElectron=
true;
521 if (!trackIsElectron) tracksToKeep.push_back(trk);
524 }
else {tracksToKeep = tracks;}
526 if(tracksToKeep.size()>=1 &&
electrons.size()>=1){
534 myPairs.push_back(pair);
551 double mass1 = massHypotheses.at(0);
552 double mass2 = massHypotheses.at(1);
562 return (mu1+mu2).M();
574 bool opposite(
false),
same(
false),
all(
false);
575 if (
selection==
"OPPOSITE") opposite=
true;
580 std::vector<JpsiEECandidate> jpsis;
583 for(
auto jpsiItr=jpsisIn.cbegin();jpsiItr!=jpsisIn.cend();jpsiItr++){
584 bool oppCh(
false),sameCh(
false);
587 qOverP2=(*jpsiItr).trackParticle2->qOverP();
588 if(qOverP1*qOverP2<0.0) oppCh=
true;
589 if(qOverP1*qOverP2>0.0) sameCh=
true;
595 tmpJpsi.
el1 = (*jpsiItr).el2;
596 tmpJpsi.
el2 = (*jpsiItr).el1;
600 if (oppCh && (opposite ||
all) ) jpsis.push_back(tmpJpsi);
601 if (sameCh && (
same ||
all) ) jpsis.push_back(tmpJpsi);
618 bool passesLHVLoose = isLHVeryLoose(*
electron);
619 bool passesLHVLoosenod0 = isLHVeryLoosenod0(*
electron);
637 return std::find(theCollection->
begin(), theCollection->
end(), theTrack) != theCollection->
end();
646 double px = 0.,
py = 0.,
pz = 0.;
647 if (0 != vxCandidate) {