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