ATLAS Offline Software
JpsiPlus2Tracks.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // ****************************************************************************
6 // ----------------------------------------------------------------------------
7 // JpsiPlus2Tracks
8 // James Catmore <James.Catmore@cern.ch>
9 // Results returned as a vector of xAOD::Vertex
10 // ----------------------------------------------------------------------------
11 // ****************************************************************************
12 
14 #include "xAODBPhys/BPhysHelper.h"
19 #include "xAODTracking/Vertex.h"
21 #include "AthLinks/ElementLink.h"
23 #include <memory>
25 #include <limits>
26 
27 namespace Analysis {
28 
29  // Set masses
30  constexpr double muMass = 105.658;
31  constexpr double kMass = 493.677;
32  constexpr double piMass = 139.57;
33  constexpr double pMass = 938.272;
34 
36 
37  // retrieving vertex Fitter
38  ATH_CHECK(m_iVertexFitter.retrieve());
39  m_VKVFitter = dynamic_cast<Trk::TrkVKalVrtFitter*>(&(*m_iVertexFitter));
40 
41  // Get the track selector tool from ToolSvc
42  ATH_CHECK(m_trkSelector.retrieve());
43 
44 
45  ATH_CHECK(m_jpsiCollectionKey.initialize());
48  if(m_MuonsUsedInJpsi.key() == "NONE") m_MuonsUsedInJpsi = "";
49  ATH_CHECK(m_MuonsUsedInJpsi.initialize(!m_MuonsUsedInJpsi.key().empty()));
50 
51  if(!m_manualMassHypo.empty() && m_manualMassHypo.size() !=4){
52  ATH_MSG_FATAL("Invalid number of elements given for manualMass hypothesis - needs 4");
53  return StatusCode::FAILURE;
54  }
55  if(!m_manualMassHypo.empty()){
56  ATH_MSG_DEBUG("manual mass hypo " << m_manualMassHypo[0] <<
57  ',' << m_manualMassHypo[1] << ',' << m_manualMassHypo[2] << ',' << m_manualMassHypo[3]);
58  }
60  ATH_MSG_FATAL("Invalid configuration");
61  return StatusCode::FAILURE;
62  }
64  ATH_MSG_FATAL("Invalid configuration");
65  return StatusCode::FAILURE;
66  }
67 
68  if(m_altMassMuonTracks.empty()){
69  m_altMassMuonTracks.assign(2, muMass); //Default to muon mass
70  }
71 
78  m_useGSFTrack.reset();
79  for(int i : m_useGSFTrackIndices) m_useGSFTrack.set(i, true);
80 
81  ATH_MSG_DEBUG("Initialize successful");
82 
83  return StatusCode::SUCCESS;
84 
85  }
86 
87  JpsiPlus2Tracks::JpsiPlus2Tracks(const std::string& t, const std::string& n, const IInterface* p) : AthAlgTool(t,n,p),
88  m_pipiMassHyp(true),
89  m_kkMassHyp(true),
90  m_kpiMassHyp(true),
91  m_kpMassHyp(false),
92  m_oppChargesOnly(true),
93  m_sameChargesOnly(false),
94  m_trkThresholdPt(0.0),
95  m_trkMaxEta(102.5),
96  m_BThresholdPt(0.0),
97  m_BMassUpper(0.0),
98  m_BMassLower(0.0),
99  m_jpsiCollectionKey("JpsiCandidates"),
100  m_jpsiMassUpper(0.0),
101  m_jpsiMassLower(0.0),
102  m_TrkParticleCollection("InDetTrackParticles"),
103  m_TrkParticleGSFCollection("GSFTrackParticles"),
104  m_MuonsUsedInJpsi(),
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),
113  m_chi2cut(-1.0),
114  m_diTrackPt(-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),
121  m_trkDeltaZ(-1.0),
122  m_requiredNMuons(0),
123  m_candidateLimit(std::numeric_limits<size_t>::max())
124  {
125  declareInterface<JpsiPlus2Tracks>(this);
126  declareProperty("pionpionHypothesis",m_pipiMassHyp);
127  declareProperty("kaonkaonHypothesis",m_kkMassHyp);
128  declareProperty("kaonpionHypothesis",m_kpiMassHyp);
129  declareProperty("kaonprotonHypothesis",m_kpMassHyp);
130  declareProperty("oppChargesOnly",m_oppChargesOnly);
131  declareProperty("SameChargesOnly",m_sameChargesOnly);
132  declareProperty("trkThresholdPt",m_trkThresholdPt);
133  declareProperty("trkMaxEta",m_trkMaxEta);
134  declareProperty("BThresholdPt",m_BThresholdPt);
135  declareProperty("BMassUpper",m_BMassUpper);
136  declareProperty("BMassLower",m_BMassLower);
137  declareProperty("JpsiContainerKey",m_jpsiCollectionKey);
138  declareProperty("JpsiMassUpper",m_jpsiMassUpper);
139  declareProperty("JpsiMassLower",m_jpsiMassLower);
140  declareProperty("TrackParticleCollection",m_TrkParticleCollection);
141  declareProperty("MuonsUsedInJpsi",m_MuonsUsedInJpsi);
142  declareProperty("ExcludeJpsiMuonsOnly",m_excludeJpsiMuonsOnly);
143  declareProperty("ExcludeCrossJpsiTracks",m_excludeCrossJpsiTracks); //Essential when trying to make vertices out of multiple muons (set to false)
144  declareProperty("TrkVertexFitterTool",m_iVertexFitter);
145  declareProperty("TrackSelectorTool",m_trkSelector);
146  declareProperty("UseMassConstraint", m_useMassConst);
147  declareProperty("AlternativeMassConstraint",m_altMassConst);
148  declareProperty("DiTrackMassUpper",m_diTrackMassUpper);
149  declareProperty("DiTrackMassLower",m_diTrackMassLower);
150  // additional cuts by Daniel Scheirich
151  declareProperty("Chi2Cut",m_chi2cut);
152  declareProperty("DiTrackPt",m_diTrackPt);
153  declareProperty("TrkQuadrupletMassUpper",m_trkQuadrupletMassUpper);
154  declareProperty("TrkQuadrupletMassLower",m_trkQuadrupletMassLower);
155  declareProperty("TrkQuadrupletPt" ,m_trkQuadrupletPt );
156  declareProperty("FinalDiTrackMassUpper" ,m_finalDiTrackMassUpper );
157  declareProperty("FinalDiTrackMassLower" ,m_finalDiTrackMassLower );
158  declareProperty("FinalDiTrackPt" ,m_finalDiTrackPt );
159  declareProperty("TrkDeltaZ" ,m_trkDeltaZ );
160  declareProperty("ManualMassHypo", m_manualMassHypo);
161  declareProperty("RequireNMuonTracks", m_requiredNMuons);
162  declareProperty("AlternativeMassConstraintTrack", m_altMassMuonTracks);
163  declareProperty("UseGSFTrackIndices", m_useGSFTrackIndices);
165  declareProperty("CandidateLimit", m_candidateLimit);
166  }
167 
169 
170 
171  //-------------------------------------------------------------------------------------
172  // Find the candidates
173  //-------------------------------------------------------------------------------------
174  StatusCode JpsiPlus2Tracks::performSearch(const EventContext& ctx, xAOD::VertexContainer& bContainer) const
175  {
176  ATH_MSG_DEBUG( "JpsiPlus2Tracks::performSearch" );
177 
178  // Get the J/psis from StoreGate
179  const xAOD::VertexContainer* importedJpsiCollection{nullptr};
181  if(!jpsihandle.isValid()){
182  ATH_MSG_ERROR("No VertexContainer with key " << m_jpsiCollectionKey.key() << " found in StoreGate. BCandidates will be EMPTY!");
183  return StatusCode::FAILURE;
184  }else{
185  importedJpsiCollection = jpsihandle.cptr();
186  ATH_MSG_DEBUG("Found VxCandidate container with key "<<m_jpsiCollectionKey.key());
187  }
188  ATH_MSG_DEBUG("VxCandidate container size " << importedJpsiCollection->size());
189 
190  // Get tracks
191  const xAOD::TrackParticleContainer* importedTrackCollection{nullptr};
193  if(!trackshandle.isValid()){
194  ATH_MSG_ERROR("No track particle collection with name " << m_TrkParticleCollection.key() << " found in StoreGate!");
195  return StatusCode::FAILURE;
196  } else {
197  importedTrackCollection = trackshandle.cptr();
198  ATH_MSG_DEBUG("Found track particle collection " << m_TrkParticleCollection.key() << " in StoreGate!");
199  }
200  ATH_MSG_DEBUG("Track container size "<< importedTrackCollection->size());
201 
202  const xAOD::TrackParticleContainer* importedGSFTrackCollection(nullptr);
203  if(m_useGSFTrack.any() && !m_TrkParticleGSFCollection.key().empty()){
205  ATH_CHECK(h.isValid());
206  importedGSFTrackCollection = h.cptr();
207  }
208 
209  // Get the muon collection used to build the J/psis
210  const xAOD::MuonContainer* importedMuonCollection{nullptr};
211  if (!m_MuonsUsedInJpsi.key().empty()) {
213  if (!h.isValid()){
214  ATH_MSG_ERROR("No muon collection with name " << m_MuonsUsedInJpsi.key() << " found in StoreGate!");
215  return StatusCode::FAILURE;
216  } else {
217  importedMuonCollection = h.cptr();
218  ATH_MSG_DEBUG("Found muon collection " << m_MuonsUsedInJpsi.key() << " in StoreGate!");
219  }
220  ATH_MSG_DEBUG("Muon container size "<< importedMuonCollection->size());
221  }
222 
223  // Typedef for vectors of tracks and VxCandidates
224  typedef std::vector<const xAOD::TrackParticle*> TrackBag;
225 
226  // Select the inner detector tracks
227  TrackBag theIDTracksAfterSelection;
228  for (auto trkPBItr=importedTrackCollection->cbegin(); trkPBItr!=importedTrackCollection->cend(); ++trkPBItr) {
229  const xAOD::TrackParticle* tp (*trkPBItr);
230  if ( tp->pt()<m_trkThresholdPt ) continue;
231  if ( std::abs(tp->eta())>m_trkMaxEta ) continue;
232  if (importedMuonCollection!=NULL && !m_excludeJpsiMuonsOnly) {
233  if (JpsiUpsilonCommon::isContainedIn(tp,importedMuonCollection)) continue;
234  }
235  if ( m_trkSelector->decision(*tp, NULL) ) theIDTracksAfterSelection.push_back(tp);
236  }
237  if (theIDTracksAfterSelection.size() == 0) return StatusCode::SUCCESS;
238  ATH_MSG_DEBUG("Number of tracks after ID trkSelector: " << theIDTracksAfterSelection.size());
239 
240  // Loop over J/psi candidates, select, collect up tracks used to build a J/psi
241  std::vector<const xAOD::Vertex*> selectedJpsiCandidates;
242  std::vector<const xAOD::TrackParticle*> jpsiTracks;
243  for(auto vxcItr=importedJpsiCollection->cbegin(); vxcItr!=importedJpsiCollection->cend(); ++vxcItr) {
244  // Check J/psi candidate invariant mass and skip if need be
245 
246  if (m_jpsiMassUpper>0.0 || m_jpsiMassLower > 0.0) {
247  xAOD::BPhysHelper jpsiCandidate(*vxcItr);
248  double jpsiMass = jpsiCandidate.totalP(m_altMassMuonTracks).M();
250  if (!pass) continue;
251  }
252  selectedJpsiCandidates.push_back(*vxcItr);
253 
254  // Collect up tracks
256  // Extract tracks from J/psi
257  const xAOD::TrackParticle* jpsiTP1 = (*vxcItr)->trackParticle(0);
258  const xAOD::TrackParticle* jpsiTP2 = (*vxcItr)->trackParticle(1);
259  jpsiTracks.push_back(jpsiTP1);
260  jpsiTracks.push_back(jpsiTP2);
261  }
262  }
263  ATH_MSG_DEBUG("selectedJpsiCandidates: " << selectedJpsiCandidates.size());
264 
265 
266 
267 
268  // Attempt to fit each track with the two tracks from the J/psi candidates
269  // Loop over J/psis
270  std::vector<const xAOD::TrackParticle*> QuadletTracks(4, nullptr);//Initialise as 4 nulls
271 
272  std::vector<double> massCuts;
273 
274  TrackBag muonTracks;
275  if (importedMuonCollection != NULL && m_excludeJpsiMuonsOnly) {
276  for(auto muon : *importedMuonCollection){
277  if(!muon->inDetTrackParticleLink().isValid()) continue;
278  auto track = muon->trackParticle( xAOD::Muon::InnerDetectorTrackParticle );
279  if(track==nullptr) continue;
280  if(!JpsiUpsilonCommon::isContainedIn(track, theIDTracksAfterSelection)) continue;
281  muonTracks.push_back(track);
282  }
283  }
284 
285  for(auto jpsiItr=selectedJpsiCandidates.begin(); jpsiItr!=selectedJpsiCandidates.end(); ++jpsiItr) {
286 
287  // Extract tracks from J/psi
288  const xAOD::TrackParticle* jpsiTP1 = (*jpsiItr)->trackParticle(0);
289  const xAOD::TrackParticle* jpsiTP2 = (*jpsiItr)->trackParticle(1);
290  QuadletTracks[0] = jpsiTP1;
291  QuadletTracks[1] = jpsiTP2;
292 
293  //If requested, only exclude duplicates in the same tripplet
295  jpsiTracks.resize(2);
296  jpsiTracks[0] = jpsiTP1;
297  jpsiTracks[1] = jpsiTP2;
298  }
299 
300  // Loop over ID tracks, call vertexing
301  for (TrackBag::iterator trkItr1=theIDTracksAfterSelection.begin(); trkItr1<theIDTracksAfterSelection.end(); ++trkItr1) { // outer loop
302  if (!m_excludeJpsiMuonsOnly && JpsiUpsilonCommon::isContainedIn(*trkItr1,jpsiTracks)) continue; // remove tracks which were used to build J/psi
303  int linkedMuonTrk1 = 0;
305  linkedMuonTrk1 = JpsiUpsilonCommon::isContainedIn(*trkItr1, muonTracks);
306  if (linkedMuonTrk1) ATH_MSG_DEBUG("This id track 1 is muon track!");
307 
308  if (JpsiUpsilonCommon::isContainedIn(*trkItr1,jpsiTracks)) {
309  if (linkedMuonTrk1) ATH_MSG_DEBUG("ID track 1 removed: id track is selected to build Jpsi!");
310  continue; // remove tracks which were used to build J/psi
311  }
312  }
313 
314  // Daniel Scheirich: remove track too far from the Jpsi vertex (DeltaZ cut)
315  if(m_trkDeltaZ>0 &&
316  std::abs((*trkItr1)->z0() + (*trkItr1)->vz() - (*jpsiItr)->z()) > m_trkDeltaZ )
317  continue;
318 
319  for (TrackBag::iterator trkItr2=trkItr1+1; trkItr2!=theIDTracksAfterSelection.end(); ++trkItr2) { // inner loop
320  if(bContainer.size() >= m_candidateLimit ){
321  ATH_MSG_WARNING("Hard Limit exceeded N=" << bContainer.size());
322  return StatusCode::SUCCESS;
323  }
324  if (!m_excludeJpsiMuonsOnly && JpsiUpsilonCommon::isContainedIn(*trkItr2,jpsiTracks)) continue; // remove tracks which were used to build J/psi
325 
327  int linkedMuonTrk2 = JpsiUpsilonCommon::isContainedIn(*trkItr2, muonTracks);
328  if (linkedMuonTrk2) ATH_MSG_DEBUG("This id track 2 is muon track!");
329  if (JpsiUpsilonCommon::isContainedIn(*trkItr2,jpsiTracks)) {
330  if (linkedMuonTrk2) ATH_MSG_DEBUG("ID track 2 removed: id track is selected to build Jpsi Vtx!");
331  continue; // remove tracks which were used to build J/psi
332  }
333  if( (linkedMuonTrk1+ linkedMuonTrk2) < m_requiredNMuons) {
334  ATH_MSG_DEBUG("Skipping Tracks with Muons " << linkedMuonTrk1 + linkedMuonTrk2 << " Limited to " << m_requiredNMuons);
335  continue;
336  }
337  }
338 
339  // Daniel Scheirich: remove track too far from the Jpsi vertex (DeltaZ cut)
340  if(m_trkDeltaZ>0 &&
341  std::abs((*trkItr2)->z0() + (*trkItr2)->vz() - (*jpsiItr)->z()) > m_trkDeltaZ )
342  continue;
343 
344  if (m_oppChargesOnly && !oppositeCharges(*trkItr1,*trkItr2)) continue; //enforce opposite charges
345  if (m_sameChargesOnly && oppositeCharges(*trkItr1,*trkItr2)) continue; //enforce same charges
346 
347  if (m_diTrackPt>0 && JpsiUpsilonCommon::getPt(*trkItr1,*trkItr2) < m_diTrackPt ) continue; // track pair pT cut (daniel Scheirich)
348 
349  // sort the track by charge, putting the positive track first
350  if((*trkItr2)->qOverP()>0) {
351  QuadletTracks[2] = *trkItr2;
352  QuadletTracks[3] = *trkItr1;
353  }else{
354  QuadletTracks[2] = *trkItr1;
355  QuadletTracks[3] = *trkItr2;
356  }
357 
358  if (m_trkQuadrupletPt>0 && JpsiUpsilonCommon::getPt(QuadletTracks[0], QuadletTracks[1], *trkItr1,*trkItr2) < m_trkQuadrupletPt ) continue; // track quadruplet pT cut (daniel Scheirich)
359 
360  // apply mass cut on track pair (non muon) if requested
361  bool passesDiTrack(true);
362  if (m_diTrackMassUpper>0.0 || m_diTrackMassLower>0.0) {
363  massCuts.clear();
364  if(m_kkMassHyp) massCuts.push_back(getInvariantMass(*trkItr1,kMass,*trkItr2,kMass));
365  if(m_pipiMassHyp) massCuts.push_back(getInvariantMass(*trkItr1,piMass,*trkItr2,piMass));
366  if(m_kpiMassHyp){
367  massCuts.push_back(getInvariantMass(*trkItr1,kMass,*trkItr2,piMass));
368  massCuts.push_back(getInvariantMass(*trkItr1,piMass,*trkItr2,kMass));
369  }
370  if(m_kpMassHyp){
371  massCuts.push_back(getInvariantMass(*trkItr1,kMass,*trkItr2,piMass));
372  massCuts.push_back(getInvariantMass(*trkItr1,piMass,*trkItr2,kMass));
373  }
374  if(!m_manualMassHypo.empty()) massCuts.push_back(getInvariantMass(*trkItr1, m_manualMassHypo[2], *trkItr2, m_manualMassHypo[3]));
376 
377  }
378  if (!passesDiTrack) continue;
379 
380  // apply mass cut on track quadruplet if requested
381  bool passes4TrackMass(true);
383  massCuts.clear();
384 
385  if(m_kkMassHyp) massCuts.push_back( getInvariantMass(QuadletTracks, m_mumukkMasses) );
386  if(m_pipiMassHyp) massCuts.push_back(getInvariantMass(QuadletTracks, m_mumupipiMasses));
387  if(m_kpiMassHyp){
388  massCuts.push_back(getInvariantMass(QuadletTracks, m_mumukpiMasses));
389  massCuts.push_back(getInvariantMass(QuadletTracks, m_mumupikMasses));
390  }
391  if(m_kpMassHyp){
392  massCuts.push_back(getInvariantMass(QuadletTracks, m_mumukpMasses));
393  massCuts.push_back(getInvariantMass(QuadletTracks, m_mumupkMasses));
394  }
395  if(!m_manualMassHypo.empty()) massCuts.push_back(getInvariantMass(QuadletTracks, m_manualMassHypo));
396 
398  }
399  if (!passes4TrackMass) continue;
400 
401  //Managed pointer, "release" if you don't want it deleted. Automatically deleted otherwise
402  std::unique_ptr<xAOD::Vertex> bVertex (fit(QuadletTracks, importedTrackCollection, importedGSFTrackCollection)); // do vertexing
403  if(!bVertex) continue;
404  double bChi2DOF = bVertex->chiSquared()/bVertex->numberDoF();
405  ATH_MSG_DEBUG("Candidate chi2/DOF is " << bChi2DOF);
406  bool chi2CutPassed = (m_chi2cut <= 0.0 || bChi2DOF < m_chi2cut);
407  if(!chi2CutPassed) { ATH_MSG_DEBUG("Chi Cut failed!"); continue; }
408  xAOD::BPhysHelper bHelper(bVertex.get());//"get" does not "release" still automatically deleted
409  bHelper.setRefTrks();
410  bool passesCuts = vertexCuts(bHelper);
411  if(!passesCuts) continue;
412  // Saving successful candidates
413  // Set links to J/psi
414  std::vector<const xAOD::Vertex*> theJpsiPreceding;
415  theJpsiPreceding.push_back(*jpsiItr);
416  bHelper.setPrecedingVertices(theJpsiPreceding, importedJpsiCollection);
417  bContainer.push_back(bVertex.release());//ptr is released preventing deletion
418 
419  } // End of inner loop over tracks
420  } // End of outer loop over tracks
421  } // End of loop over J/spi
422  ATH_MSG_DEBUG("bContainer size " << bContainer.size());
423  return StatusCode::SUCCESS;
424 
425  }
426 
427  // *********************************************************************************
428 
429  // ---------------------------------------------------------------------------------
430  // fit - does the fit
431  // ---------------------------------------------------------------------------------
432 
433  xAOD::Vertex* JpsiPlus2Tracks::fit(const std::vector<const xAOD::TrackParticle*> &inputTracks,
434  const xAOD::TrackParticleContainer* importedTrackCollection,
435  const xAOD::TrackParticleContainer* gsfCollection) const {
436 
437  std::unique_ptr<Trk::IVKalState> state = m_VKVFitter->makeState();
438 
439 
440 
441 
442  // Set the mass constraint if requested by user (default=true)
443  // Can be set by user (m_altMassConstraint) - default is -1.0.
444  // If < 0.0, uses J/psi (default)
445  // If > 0.0, uses the value provided
446 
447 
448  if (m_useMassConst) {
449  constexpr double jpsiTableMass = 3096.916;
451  std::array<int,2> indices= {1, 2};
452  if (m_altMassConst<0.0) m_VKVFitter->setMassForConstraint(jpsiTableMass,indices,*state);
454  }
455 
456  // Do the fit itself.......
457  // Starting point (use the J/psi position)
458  Amg::Vector3D startingPoint(0,0,0);
459  StatusCode sc=m_VKVFitter->VKalVrtFitFast(inputTracks, startingPoint, *state);
460  if(sc.isFailure()){
461  startingPoint = Amg::Vector3D(0,0,0);
462  }
463  xAOD::Vertex* theResult = m_VKVFitter->fit(inputTracks, startingPoint, *state);
464 
465  // Added by ASC
466  if(theResult != 0){
467  std::vector<ElementLink<DataVector<xAOD::TrackParticle> > > newLinkVector;
468  for(unsigned int i=0; i< theResult->trackParticleLinks().size(); i++)
469  {
470  ElementLink<DataVector<xAOD::TrackParticle> > mylink=theResult->trackParticleLinks()[i]; //makes a copy (non-const)
471  mylink.setStorableObject( m_useGSFTrack[i] ? *gsfCollection : *importedTrackCollection, true);
472  newLinkVector.push_back( mylink );
473  }
474  theResult->clearTracks();
475  theResult->setTrackParticleLinks( newLinkVector );
476  }
477 
478  return theResult;
479 
480  }
481 
482 
483 
484  // *********************************************************************************
485 
486  // oppositeCharges: true if the two tracks are oppositely charged
488  double qOverP1=trk1->qOverP();
489  double qOverP2=trk2->qOverP();
490  bool opposite = (qOverP1*qOverP2<0.0);
491  return opposite;
492  }
493 
494 
495 
496  // ---------------------------------------------------------------------------------
497  // getInvariantMass: returns invariant mass given a pair of tracks and their mass
498  // hypothesis. Each track must have a separate mass hypothesis in
499  // the vector, and they must be in the same order as the tracks in the track vector.
500  // Otherwise it will go horribly wrong.
501  // ---------------------------------------------------------------------------------
502 
503  double JpsiPlus2Tracks::getInvariantMass(const xAOD::TrackParticle* trk1, double mass1, const xAOD::TrackParticle* trk2, double mass2){
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;
517  double eSum=e1+e2;
518  double M=sqrt((eSum*eSum)-(pxSum*pxSum)-(pySum*pySum)-(pzSum*pzSum));
519 
520  return M;
521 
522  }
523 
524  double JpsiPlus2Tracks::getInvariantMass(const std::vector<const xAOD::TrackParticle*> &trk,
525  const std::vector<double> &masses)
526  {
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();
532  double e1 = sqrt(px1*px1+py1*py1+pz1*pz1+masses[0]*masses[0]);
533 
534  const auto trk2V = trk[1]->p4();
535  double px2 = trk2V.Px();
536  double py2 = trk2V.Py();
537  double pz2 = trk2V.Pz();
538  double e2 = sqrt(px2*px2+py2*py2+pz2*pz2+masses[1]*masses[1]);
539 
540  const auto trk3V = trk[2]->p4();
541  double px3 = trk3V.Px();
542  double py3 = trk3V.Py();
543  double pz3 = trk3V.Pz();
544  double e3 = sqrt(px3*px3+py3*py3+pz3*pz3+masses[2]*masses[2]);
545 
546  const auto trk4V = trk[3]->p4();
547  double px4 = trk4V.Px();
548  double py4 = trk4V.Py();
549  double pz4 = trk4V.Pz();
550  double e4 = sqrt(px4*px4+py4*py4+pz4*pz4+masses[3]*masses[3]);
551 
552  double pxSum=px1+px2+px3+px4;
553  double pySum=py1+py2+py3+py4;
554  double pzSum=pz1+pz2+pz3+pz4;
555  double eSum=e1+e2+e3+e4;
556 
557  double M=sqrt((eSum*eSum)-(pxSum*pxSum)-(pySum*pySum)-(pzSum*pzSum));
558 
559  return M;
560 
561  }
562 
563  bool JpsiPlus2Tracks::passCuts(xAOD::BPhysHelper &bHelper, std::span<const double> masses, std::string_view str) const{
564 
565  TLorentzVector bMomentum = bHelper.totalP(masses);
566 // ATH_MSG_DEBUG(bMomentum.X() << " " << bMomentum.Y()<< " " << bMomentum.Z() << " " << bMomentum.E());
567  double bPt = bMomentum.Pt();
568 
569  double bMass = bMomentum.M();
570  ATH_MSG_DEBUG("Candidate pt/mass under " << str << " track mass hypothesis is " << bPt << " / " << bMass);
571  if( !JpsiUpsilonCommon::cutAcceptGreater(bPt, m_BThresholdPt)) return false;
572  if( !JpsiUpsilonCommon::cutRange(bMass, m_BMassLower, m_BMassUpper)) return false;
573  assert(masses.size()==4);
574  TLorentzVector tr1 = bHelper.refTrk(2,masses[2]);
575  TLorentzVector tr2 = bHelper.refTrk(3,masses[3]);
576  double bDiTrkPt = (tr1+tr2).Pt();
577  double bDiTrkMass = (tr1+tr2).M();
578  if( !JpsiUpsilonCommon::cutAcceptGreater(bDiTrkPt, m_finalDiTrackPt)) return false;
580 
581  return true;
582  }
583 
585  // Invariant mass calculation under the various hypotheses
586  // create the helper class
587 
588  bool passesCuts = (m_kkMassHyp && passCuts(bHelper, m_mumukkMasses, "KK")) ||
589  (m_pipiMassHyp && passCuts(bHelper, m_mumupipiMasses, "pi pi")) ||
590  (m_kpiMassHyp && (passCuts(bHelper, m_mumukpiMasses, "K pi") ||
591  passCuts(bHelper, m_mumupikMasses, "pi K"))) ||
592  (m_kpMassHyp && (passCuts(bHelper, m_mumukpMasses, "K p") ||
593  passCuts(bHelper, m_mumupkMasses, "p K"))) ||
594  (!m_manualMassHypo.empty() && passCuts(bHelper, m_manualMassHypo, "manual"));
595  return passesCuts;
596  }
597 
598 
599 } // End of namespace
600 
601 
602 
Analysis::JpsiPlus2Tracks::m_oppChargesOnly
bool m_oppChargesOnly
Definition: JpsiPlus2Tracks.h:71
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
Analysis::JpsiPlus2Tracks::m_jpsiMassLower
double m_jpsiMassLower
Definition: JpsiPlus2Tracks.h:80
Analysis::JpsiPlus2Tracks::m_trkSelector
ToolHandle< Trk::ITrackSelectorTool > m_trkSelector
Definition: JpsiPlus2Tracks.h:87
Analysis::JpsiPlus2Tracks::m_MuonsUsedInJpsi
SG::ReadHandleKey< xAOD::MuonContainer > m_MuonsUsedInJpsi
Definition: JpsiPlus2Tracks.h:83
xAOD::muon
@ muon
Definition: TrackingPrimitives.h:195
xAOD::BPhysHelper::totalP
TVector3 totalP()
: Returns total 3-momentum calculated from the refitted tracks
Definition: BPhysHelper.cxx:374
Analysis::JpsiPlus2Tracks::m_BMassLower
double m_BMassLower
Definition: JpsiPlus2Tracks.h:77
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
Analysis::JpsiPlus2Tracks::m_chi2cut
double m_chi2cut
Definition: JpsiPlus2Tracks.h:95
Analysis::JpsiPlus2Tracks::~JpsiPlus2Tracks
~JpsiPlus2Tracks()
Definition: JpsiPlus2Tracks.cxx:168
Analysis::JpsiPlus2Tracks::m_kkMassHyp
bool m_kkMassHyp
Definition: JpsiPlus2Tracks.h:68
Analysis::JpsiPlus2Tracks::JpsiPlus2Tracks
JpsiPlus2Tracks(const std::string &t, const std::string &n, const IInterface *p)
Definition: JpsiPlus2Tracks.cxx:87
xAOD::BPhysHelper
Definition: BPhysHelper.h:71
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
max
#define max(a, b)
Definition: cfImp.cxx:41
Analysis::JpsiPlus2Tracks::m_BThresholdPt
double m_BThresholdPt
Definition: JpsiPlus2Tracks.h:75
VertexPointEstimator.h
SG::ReadHandle::cptr
const_pointer_type cptr()
Dereference the pointer.
Analysis::JpsiPlus2Tracks::m_trkQuadrupletPt
double m_trkQuadrupletPt
Definition: JpsiPlus2Tracks.h:99
egammaEnergyPositionAllSamples::e1
double e1(const xAOD::CaloCluster &cluster)
return the uncorrected cluster energy in 1st sampling
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:70
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
Analysis::JpsiPlus2Tracks::m_finalDiTrackMassLower
double m_finalDiTrackMassLower
Definition: JpsiPlus2Tracks.h:101
Trk::TrkVKalVrtFitter::VKalVrtFitFast
virtual StatusCode VKalVrtFitFast(std::span< const xAOD::TrackParticle *const >, Amg::Vector3D &Vertex, double &minDZ, IVKalState &istate) const
Definition: VKalVrtFitFastSvc.cxx:56
xAOD::Vertex_v1::trackParticleLinks
const TrackParticleLinks_t & trackParticleLinks() const
Get all the particles associated with the vertex.
Analysis::JpsiPlus2Tracks::m_candidateLimit
size_t m_candidateLimit
Definition: JpsiPlus2Tracks.h:116
Trk::indices
std::pair< long int, long int > indices
Definition: AlSymMatBase.h:24
Analysis::JpsiPlus2Tracks::m_kpiMassHyp
bool m_kpiMassHyp
Definition: JpsiPlus2Tracks.h:69
ParticleTest.tp
tp
Definition: ParticleTest.py:25
Analysis::kMass
constexpr double kMass
Definition: JpsiPlus1Track.cxx:34
Analysis::JpsiPlus2Tracks::m_mumupikMasses
std::vector< double > m_mumupikMasses
Definition: JpsiPlus2Tracks.h:111
Analysis::JpsiPlus2Tracks::fit
xAOD::Vertex * fit(const std::vector< const xAOD::TrackParticle * > &, const xAOD::TrackParticleContainer *, const xAOD::TrackParticleContainer *GSL) const
Definition: JpsiPlus2Tracks.cxx:433
Analysis::JpsiPlus2Tracks::initialize
virtual StatusCode initialize() override
Definition: JpsiPlus2Tracks.cxx:35
Analysis::JpsiUpsilonCommon::isContainedIn
static bool isContainedIn(const xAOD::TrackParticle *, const std::vector< const xAOD::TrackParticle * > &)
Definition: JpsiUpsilonCommon.cxx:58
xAOD::BPhysHelper::setPrecedingVertices
bool setPrecedingVertices(const std::vector< const xAOD::Vertex * > &vertices, const xAOD::VertexContainer *vertexContainer)
Sets links to preceding vertices.
Definition: BPhysHelper.cxx:641
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
SG::VarHandleKey::key
const std::string & key() const
Return the StoreGate ID for the referenced object.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:141
Analysis::pMass
constexpr double pMass
Definition: JpsiPlus2Tracks.cxx:33
Analysis::JpsiPlus2Tracks::m_sameChargesOnly
bool m_sameChargesOnly
Definition: JpsiPlus2Tracks.h:72
Analysis::JpsiPlus2Tracks::m_trkQuadrupletMassLower
double m_trkQuadrupletMassLower
Definition: JpsiPlus2Tracks.h:98
Analysis::JpsiPlus2Tracks::m_VKVFitter
Trk::TrkVKalVrtFitter * m_VKVFitter
Definition: JpsiPlus2Tracks.h:88
Analysis::JpsiPlus2Tracks::performSearch
virtual StatusCode performSearch(const EventContext &ctx, xAOD::VertexContainer &) const override
Definition: JpsiPlus2Tracks.cxx:174
Analysis::JpsiUpsilonCommon::cutRange
static bool cutRange(double value, double min, double max)
Definition: JpsiUpsilonCommon.cxx:72
Analysis::JpsiPlus2Tracks::oppositeCharges
static bool oppositeCharges(const xAOD::TrackParticle *, const xAOD::TrackParticle *)
Definition: JpsiPlus2Tracks.cxx:487
TrkVKalVrtFitter.h
AthenaPoolTestRead.sc
sc
Definition: AthenaPoolTestRead.py:27
Analysis::JpsiPlus2Tracks::m_finalDiTrackPt
double m_finalDiTrackPt
Definition: JpsiPlus2Tracks.h:102
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
Analysis::JpsiPlus2Tracks::m_mumukkMasses
std::vector< double > m_mumukkMasses
Definition: JpsiPlus2Tracks.h:108
xAOD::TrackParticle_v1::p4
virtual FourMom_t p4() const override final
The full 4-momentum of the particle.
Definition: TrackParticle_v1.cxx:129
Analysis::JpsiPlus2Tracks::m_mumupipiMasses
std::vector< double > m_mumupipiMasses
Definition: JpsiPlus2Tracks.h:109
Analysis::JpsiPlus2Tracks::m_TrkParticleCollection
SG::ReadHandleKey< xAOD::TrackParticleContainer > m_TrkParticleCollection
Definition: JpsiPlus2Tracks.h:81
Analysis::JpsiUpsilonCommon::cutAcceptGreater
static bool cutAcceptGreater(double value, double min)
Definition: JpsiUpsilonCommon.cxx:83
Trk::TrkVKalVrtFitter::setMassForConstraint
virtual void setMassForConstraint(double Mass, IVKalState &istate) const override final
Definition: SetFitOptions.cxx:134
Analysis::JpsiPlus2Tracks::getInvariantMass
static double getInvariantMass(const xAOD::TrackParticle *, double, const xAOD::TrackParticle *, double)
Definition: JpsiPlus2Tracks.cxx:503
Analysis::JpsiPlus2Tracks::m_pipiMassHyp
bool m_pipiMassHyp
Definition: JpsiPlus2Tracks.h:67
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
CheckAppliedSFs.e3
e3
Definition: CheckAppliedSFs.py:264
python.TrigInDetConfig.inputTracks
inputTracks
Definition: TrigInDetConfig.py:168
lumiFormat.i
int i
Definition: lumiFormat.py:92
PixelAthClusterMonAlgCfg.e4
e4
Definition: PixelAthClusterMonAlgCfg.py:317
beamspotman.n
n
Definition: beamspotman.py:731
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
extractSporadic.h
list h
Definition: extractSporadic.py:97
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
Analysis::JpsiPlus2Tracks::m_mumukpMasses
std::vector< double > m_mumukpMasses
Definition: JpsiPlus2Tracks.h:112
Analysis::JpsiPlus2Tracks::m_jpsiMassUpper
double m_jpsiMassUpper
Definition: JpsiPlus2Tracks.h:79
xAOD::Vertex_v1::setTrackParticleLinks
void setTrackParticleLinks(const TrackParticleLinks_t &trackParticles)
Set all track particle links at once.
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Analysis::JpsiPlus2Tracks::m_altMassConst
double m_altMassConst
Definition: JpsiPlus2Tracks.h:90
Analysis::JpsiPlus2Tracks::m_excludeCrossJpsiTracks
bool m_excludeCrossJpsiTracks
Definition: JpsiPlus2Tracks.h:85
JpsiUpsilonCommon.h
Analysis::JpsiPlus2Tracks::m_diTrackPt
double m_diTrackPt
Definition: JpsiPlus2Tracks.h:96
DerivationFramework::TrackBag
std::vector< const xAOD::TrackParticle * > TrackBag
Definition: BPhysAddMuonBasedInvMass.h:32
xAOD::Vertex_v1::clearTracks
void clearTracks()
Remove all tracks from the vertex.
Definition: Vertex_v1.cxx:331
SG::VarHandleKey::initialize
StatusCode initialize(bool used=true)
If this object is used as a property, then this should be called during the initialize phase.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:103
Analysis::JpsiPlus2Tracks::m_finalDiTrackMassUpper
double m_finalDiTrackMassUpper
Definition: JpsiPlus2Tracks.h:100
DataVector
Derived DataVector<T>.
Definition: DataVector.h:581
Vertex.h
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
Analysis::JpsiPlus2Tracks::m_useGSFTrack
std::bitset< 4 > m_useGSFTrack
Definition: JpsiPlus2Tracks.h:115
Analysis::JpsiPlus2Tracks::m_excludeJpsiMuonsOnly
bool m_excludeJpsiMuonsOnly
Definition: JpsiPlus2Tracks.h:84
Analysis::JpsiPlus2Tracks::m_altMassMuonTracks
std::vector< double > m_altMassMuonTracks
Definition: JpsiPlus2Tracks.h:107
Analysis::JpsiPlus2Tracks::vertexCuts
bool vertexCuts(xAOD::BPhysHelper &bHelper) const
Definition: JpsiPlus2Tracks.cxx:584
xAOD::TrackParticle_v1::qOverP
float qOverP() const
Returns the parameter.
Analysis
The namespace of all packages in PhysicsAnalysis/JetTagging.
Definition: BTaggingCnvAlg.h:20
JpsiPlus2Tracks.h
Analysis::JpsiPlus2Tracks::m_BMassUpper
double m_BMassUpper
Definition: JpsiPlus2Tracks.h:76
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
Analysis::JpsiPlus2Tracks::m_jpsiCollectionKey
SG::ReadHandleKey< xAOD::VertexContainer > m_jpsiCollectionKey
Definition: JpsiPlus2Tracks.h:78
Analysis::JpsiPlus2Tracks::m_requiredNMuons
int m_requiredNMuons
Definition: JpsiPlus2Tracks.h:106
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Analysis::JpsiPlus2Tracks::m_TrkParticleGSFCollection
SG::ReadHandleKey< xAOD::TrackParticleContainer > m_TrkParticleGSFCollection
Definition: JpsiPlus2Tracks.h:82
IVertexFitter.h
Analysis::JpsiPlus2Tracks::m_useGSFTrackIndices
std::vector< int > m_useGSFTrackIndices
Definition: JpsiPlus2Tracks.h:114
BPhysHelper.h
: B-physics xAOD helpers.
xAOD::Vertex_v1::numberDoF
float numberDoF() const
Returns the number of degrees of freedom of the vertex fit as float.
TrackParticle.h
xAOD::BPhysHelper::setRefTrks
bool setRefTrks(std::vector< float > px, std::vector< float > py, std::vector< float > pz)
Sets refitted track momenta.
Definition: BPhysHelper.cxx:286
VertexContainer.h
Prompt::Def::Pt
@ Pt
Definition: VarHolder.h:76
Analysis::JpsiPlus2Tracks::m_useMassConst
bool m_useMassConst
Definition: JpsiPlus2Tracks.h:89
xAOD::Vertex_v1::chiSquared
float chiSquared() const
Returns the of the vertex fit as float.
Analysis::JpsiPlus2Tracks::m_mumupkMasses
std::vector< double > m_mumupkMasses
Definition: JpsiPlus2Tracks.h:113
h
xAOD::Vertex_v1
Class describing a Vertex.
Definition: Vertex_v1.h:42
xAOD::BPhysHelper::refTrk
TVector3 refTrk(const size_t index)
Returns i-th refitted track 3-momentum.
Definition: BPhysHelper.cxx:126
Analysis::JpsiPlus2Tracks::m_trkMaxEta
double m_trkMaxEta
Definition: JpsiPlus2Tracks.h:74
egammaEnergyPositionAllSamples::e2
double e2(const xAOD::CaloCluster &cluster)
return the uncorrected cluster energy in 2nd sampling
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Analysis::JpsiPlus2Tracks::m_kpMassHyp
bool m_kpMassHyp
Definition: JpsiPlus2Tracks.h:70
Analysis::piMass
constexpr double piMass
Definition: JpsiPlus1Track.cxx:35
Analysis::JpsiPlus2Tracks::passCuts
bool passCuts(xAOD::BPhysHelper &bHelper, std::span< const double > masses, std::string_view str) const
Definition: JpsiPlus2Tracks.cxx:563
Analysis::JpsiPlus2Tracks::m_trkQuadrupletMassUpper
double m_trkQuadrupletMassUpper
Definition: JpsiPlus2Tracks.h:97
str
Definition: BTagTrackIpAccessor.cxx:11
xAOD::track
@ track
Definition: TrackingPrimitives.h:512
xAOD::TrackParticle_v1
Class describing a TrackParticle.
Definition: TrackParticle_v1.h:43
Trk::TrkVKalVrtFitter::fit
virtual xAOD::Vertex * fit(const std::vector< const TrackParameters * > &perigeeList, const Amg::Vector3D &startingPoint) const override final
Interface for MeasuredPerigee with starting point.
Definition: TrkVKalVrtFitter.cxx:261
Analysis::JpsiPlus2Tracks::m_iVertexFitter
ToolHandle< Trk::IVertexFitter > m_iVertexFitter
Definition: JpsiPlus2Tracks.h:86
AthAlgTool
Definition: AthAlgTool.h:26
Analysis::muMass
constexpr double muMass
Definition: JpsiPlus1Track.cxx:33
Analysis::JpsiPlus2Tracks::m_trkDeltaZ
double m_trkDeltaZ
Definition: JpsiPlus2Tracks.h:103
ITrackSelectorTool.h
Analysis::JpsiPlus2Tracks::m_trkThresholdPt
double m_trkThresholdPt
Definition: JpsiPlus2Tracks.h:73
Analysis::JpsiUpsilonCommon::getPt
static double getPt(const xAOD::TrackParticle *, const xAOD::TrackParticle *)
Definition: JpsiUpsilonCommon.cxx:17
Analysis::JpsiPlus2Tracks::m_mumukpiMasses
std::vector< double > m_mumukpiMasses
Definition: JpsiPlus2Tracks.h:110
Analysis::JpsiPlus2Tracks::m_manualMassHypo
std::vector< double > m_manualMassHypo
Definition: JpsiPlus2Tracks.h:105
DataVector::size
size_type size() const noexcept
Returns the number of elements in the collection.
Trk::TrkVKalVrtFitter::setMassInputParticles
virtual void setMassInputParticles(const std::vector< double > &, IVKalState &istate) const override final
Definition: SetFitOptions.cxx:187
Analysis::JpsiPlus2Tracks::m_diTrackMassUpper
double m_diTrackMassUpper
Definition: JpsiPlus2Tracks.h:91
Trk::TrkVKalVrtFitter::makeState
virtual std::unique_ptr< IVKalState > makeState(const EventContext &ctx) const override final
Definition: TrkVKalVrtFitter.cxx:118
SUSY_SimplifiedModel_PreInclude.masses
dictionary masses
Definition: SUSY_SimplifiedModel_PreInclude.py:7
Analysis::JpsiPlus2Tracks::m_diTrackMassLower
double m_diTrackMassLower
Definition: JpsiPlus2Tracks.h:92
Analysis::JpsiUpsilonCommon::cutRangeOR
static bool cutRangeOR(const std::vector< double > &values, double min, double max)
Definition: JpsiUpsilonCommon.cxx:76
Trk::TrkVKalVrtFitter
Definition: TrkVKalVrtFitter.h:67