ATLAS Offline Software
MuPlusDsCascade.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
5 // MuPlusDsCascade.cxx, (c) ATLAS Detector software
15 
16 #include <algorithm>
17 #include "HepPDT/ParticleDataTable.hh"
18 
19 #include "xAODMuon/MuonContainer.h"
24 
29 
30 namespace DerivationFramework {
32  typedef std::vector<VertexLink> VertexLinkVector;
33  typedef std::vector<const xAOD::TrackParticle*> TrackBag;
34 
36  typedef std::vector<MuonsLink> MuonsLinkVector;
37 
39 
40  // retrieving vertex Fitter
41  ATH_CHECK( m_iVertexFitter.retrieve());
42 
43  // retrieving the V0 tools //from JpsiFinder
44  ATH_CHECK( m_V0Tools.retrieve());
45 
46  // Get the track selector tool from ToolSvc
47  if ( m_trkSelector.retrieve().isFailure() ) {
48  ATH_MSG_FATAL("Failed to retrieve tool " << m_trkSelector);
49  return StatusCode::FAILURE;
50  } else {
51  ATH_MSG_INFO("Retrieved tool " << m_trkSelector);
52  }
53 
54  //======================== inDetTrack selection tool ==================
55  m_trackSelectionTools = std::make_unique<InDet::InDetTrackSelectionTool>("TrackSelector");
56  ANA_CHECK(m_trackSelectionTools->setProperty("CutLevel", "LoosePrimary"));
57  ANA_CHECK(m_trackSelectionTools->initialize() );
58  //=====================================================================
59 
60  // retrieving the Cascade tools
61  ATH_CHECK( m_CascadeTools.retrieve());
62 
64 
65  ATH_CHECK( m_partPropSvc.retrieve() );
66  auto pdt = m_partPropSvc->PDT();
67 
68  // Ds+/- : K K π
69  // D+/- : K π π
70  // Lambda_c+/- : K π p
71 
72  // retrieve particle masses
73  if(m_vtx0MassHypo < 0.)
75  if(m_vtx1MassHypo < 0.) {
76  if(std::abs(m_Dx_pid) == 411) m_vtx1MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::DPLUS);
77  if(std::abs(m_Dx_pid) == 431) m_vtx1MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::DSPLUS);
78  if(std::abs(m_Dx_pid) == 4122) m_vtx1MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::LAMBDACPLUS);
79  }
80 
82  if(m_vtx1Daug1MassHypo < 0.) {
83  if(std::abs(m_Dx_pid) == 431) m_vtx1Daug1MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::KPLUS); //Ds+
84  else m_vtx1Daug1MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::PIPLUS); //D+, Lambda_c+
85  }
87  if(m_vtx1Daug3MassHypo < 0.) {
88  if(std::abs(m_Dx_pid) == 4122) m_vtx1Daug3MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::PROTON); //Lambda_c+
89  else m_vtx1Daug3MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::PIPLUS); //Ds+, D+
90  }
91  return StatusCode::SUCCESS;
92  }
93 
94 
96  {
97  std::vector<Trk::VxCascadeInfo*> cascadeinfoContainer;
98  constexpr int topoN = 2;
99  std::array<xAOD::VertexContainer*, topoN> Vtxwritehandles;
100  std::array<xAOD::VertexAuxContainer*, topoN> Vtxwritehandlesaux;
101  if(m_cascadeOutputsKeys.size() !=topoN) { ATH_MSG_FATAL("Incorrect number of VtxContainers"); return StatusCode::FAILURE; }
102 
103  for(int i =0; i<topoN;i++){
104  Vtxwritehandles[i] = new xAOD::VertexContainer();
105  Vtxwritehandlesaux[i] = new xAOD::VertexAuxContainer();
106  Vtxwritehandles[i]->setStore(Vtxwritehandlesaux[i]);
107  ATH_CHECK(evtStore()->record(Vtxwritehandles[i] , m_cascadeOutputsKeys[i] ));
108  ATH_CHECK(evtStore()->record(Vtxwritehandlesaux[i], m_cascadeOutputsKeys[i] + "Aux."));
109  }
110 
111  //----------------------------------------------------
112  // retrieve primary vertices
113  //----------------------------------------------------
114  const xAOD::Vertex * primaryVertex(nullptr);
115  const xAOD::VertexContainer *pvContainer(nullptr);
117  ATH_MSG_DEBUG("Found " << m_VxPrimaryCandidateName << " in StoreGate!");
118 
119  if (pvContainer->size()==0){
120  ATH_MSG_WARNING("You have no primary vertices: " << pvContainer->size());
121  return StatusCode::RECOVERABLE;
122  } else {
123  primaryVertex = (*pvContainer)[0];
124  }
125 
126  //----------------------------------------------------
127  // Try to retrieve refitted primary vertices
128  //----------------------------------------------------
129  xAOD::VertexContainer* refPvContainer = nullptr;
130  xAOD::VertexAuxContainer* refPvAuxContainer = nullptr;
131  if (m_refitPV) {
132  if (evtStore()->contains<xAOD::VertexContainer>(m_refPVContainerName)) {
133  // refitted PV container exists. Get it from the store gate
134  ATH_CHECK(evtStore()->retrieve(refPvContainer , m_refPVContainerName ));
135  ATH_CHECK(evtStore()->retrieve(refPvAuxContainer, m_refPVContainerName + "Aux."));
136  } else {
137  // refitted PV container does not exist. Create a new one.
138  refPvContainer = new xAOD::VertexContainer;
139  refPvAuxContainer = new xAOD::VertexAuxContainer;
140  refPvContainer->setStore(refPvAuxContainer);
141  ATH_CHECK(evtStore()->record(refPvContainer , m_refPVContainerName));
142  ATH_CHECK(evtStore()->record(refPvAuxContainer, m_refPVContainerName+"Aux."));
143  }
144  }
145 
146  ATH_CHECK(performSearch(&cascadeinfoContainer));
148  if(!evt.isValid()) {
149  ATH_MSG_ERROR("Cannot Retrieve " << m_eventInfo_key.key() );
150  return StatusCode::FAILURE;
151  }
152 
154  helper.SetMinNTracksInPV(m_PV_minNTracks);
155 
156  // Decorators for the main vertex: chi2, ndf, pt and pt error, plus the V0 vertex variables
157  SG::AuxElement::Decorator<VertexLinkVector> CascadeLinksDecor("CascadeVertexLinks");
158  SG::AuxElement::Decorator<MuonsLinkVector> MuonsLinksDecor("MuonsLinks");
159  SG::AuxElement::Decorator<VertexLinkVector> DxLinksDecor("DxVertexLinks");
160  SG::AuxElement::Decorator<float> chi2_decor("ChiSquared");
161  SG::AuxElement::Decorator<float> ndof_decor("NumberDoF");
162  SG::AuxElement::Decorator<float> Pt_decor("Pt");
163  SG::AuxElement::Decorator<float> PtErr_decor("PtErr");
164  SG::AuxElement::Decorator<float> Mass_svdecor("Dx_mass");
165  SG::AuxElement::Decorator<float> MassErr_svdecor("Dx_massErr");
166  SG::AuxElement::Decorator<float> Pt_svdecor("Dx_Pt");
167  SG::AuxElement::Decorator<float> PtErr_svdecor("Dx_PtErr");
168  SG::AuxElement::Decorator<float> Lxy_svdecor("Dx_Lxy");
169  SG::AuxElement::Decorator<float> LxyErr_svdecor("Dx_LxyErr");
170  SG::AuxElement::Decorator<float> Tau_svdecor("Dx_Tau");
171  SG::AuxElement::Decorator<float> TauErr_svdecor("Dx_TauErr");
172  SG::AuxElement::Decorator<float> Dx_chi2_svdecor("Dx_chi2");
173  SG::AuxElement::Decorator<float> Dx_ndof_svdecor("Dx_ndof");
174 
175  // K X1 X2
176  // Ds+/- : K K π
177  // D+/- : K π π
178  // Lambda_c+/- : K π p
179 
180  SG::AuxElement::Decorator<float> massKX1_svdecor("KX_mass");
181  SG::AuxElement::Decorator<float> massKX1X2_svdecor("KXpi_mass");
182  SG::AuxElement::Decorator<float> RapidityKX1X2_svdecor("KXpi_Rapidity");
183 
184  SG::AuxElement::Decorator<float> MuMass_decor("Mu_mass");
185  SG::AuxElement::Decorator<float> MuPt_decor("Mu_pt");
186  SG::AuxElement::Decorator<float> MuEta_decor("Mu_eta");
187  SG::AuxElement::Decorator<float> MuChi2_decor("Mu_chi2");
188  SG::AuxElement::Decorator<float> MunDoF_decor("Mu_nDoF");
189  //muon contribution to chi2 of the cascade fit
190  SG::AuxElement::Decorator<float> MuChi2B_decor("Mu_chi2_B");
191  SG::AuxElement::Decorator<float> MunDoFB_decor("Mu_nDoF_B");
192 
193  SG::AuxElement::Decorator<float> ChargeK_decor("K_charge");
194  SG::AuxElement::Decorator<float> ChargeX1_decor("X_charge"); // K (Ds+) or pi (D+, Lambda_c+)
195  SG::AuxElement::Decorator<float> ChargeX2_decor("Pi_charge"); // pi (Ds+, D+) or p (Lambda_c+)
196  SG::AuxElement::Decorator<float> ChargeMu_decor("Mu_charge");
197 
198  ATH_MSG_DEBUG("cascadeinfoContainer size " << cascadeinfoContainer.size());
199 
200  // Get Muons container
201  const xAOD::MuonContainer* muonContainer(nullptr);
203  ATH_MSG_DEBUG("Muon container size "<<muonContainer->size());
204 
205  // Get D_(s)+/Lambda_c+ container and identify the input D_(s)+/Lambda_c+
206  const xAOD::VertexContainer *dxContainer(nullptr);
208 
209  for (Trk::VxCascadeInfo* x : cascadeinfoContainer) {
210  if(x==nullptr){
211  ATH_MSG_ERROR("cascadeinfoContainer is null");
212  continue;
213  }
214 
215  // the cascade fitter returns:
216  // std::vector<xAOD::Vertex*>, each xAOD::Vertex contains the refitted track parameters (perigee at the vertex position)
217  // vertices[iv] the links to the original TPs and a covariance of size 3+5*NTRK; the chi2 of the total fit
218  // is split between the cascade vertices as per track contribution
219  // std::vector< std::vector<TLorentzVector> >, each std::vector<TLorentzVector> contains the refitted momenta (TLorentzVector)
220  // momenta[iv][...] of all tracks in the corresponding vertex, including any pseudotracks (from cascade vertices)
221  // originating in this vertex; the masses are as assigned in the cascade fit
222  // std::vector<Amg::MatrixX>, the corresponding covariance matrices in momentum space
223  // covariance[iv]
224  // int nDoF, double Chi2
225  //
226  // the invariant mass, pt, lifetime etc. errors should be calculated using the covariance matrices in momentum space as these
227  // take into account the full track-track and track-vertex correlations
228  //
229  // in the case of Jpsi+V0: vertices[0] is the V0 vertex, vertices[1] is the B/Lambda_b(bar) vertex, containing the 2 Jpsi tracks.
230  // The covariance terms between the two vertices are not stored. In momentum space momenta[0] contains the 2 V0 tracks,
231  // their momenta add up to the momentum of the 3rd track in momenta[1], the first two being the Jpsi tracks
232 
233  const std::vector<xAOD::Vertex*> &cascadeVertices = x->vertices();
234 
235  if(cascadeVertices.size()!=topoN) ATH_MSG_ERROR("Incorrect number of vertices");
236  if(cascadeVertices[0] == nullptr || cascadeVertices[1] == nullptr) ATH_MSG_ERROR("Error null vertex");
237  // Keep vertices (bear in mind that they come in reverse order!)
238  for(int i =0;i<topoN;i++) Vtxwritehandles[i]->push_back(cascadeVertices[i]);
239 
240  x->setSVOwnership(false); // Prevent Container from deleting vertices
241  const auto mainVertex = cascadeVertices[1]; // this is the B_c+/- vertex
242  //and cascadeVertices[0] is Dx vertex
243  const std::vector< std::vector<TLorentzVector> > &moms = x->getParticleMoms();
244 
245 
246  // Set links to cascade vertices
247  std::vector<const xAOD::Vertex*> verticestoLink;
248  verticestoLink.push_back(cascadeVertices[0]);
249  if(Vtxwritehandles[1] == nullptr) ATH_MSG_ERROR("Vtxwritehandles[1] is null");
250  if(!BPhysPVCascadeTools::LinkVertices(CascadeLinksDecor, verticestoLink, Vtxwritehandles[0], cascadeVertices[1]))
251  ATH_MSG_ERROR("Error decorating with cascade vertices");
252 
253  //----------------------------------- mu
254  // Identify muons
255  typedef std::vector<const xAOD::Muon*> MuonBag;
256  MuonBag selectedMuons; selectedMuons.clear();
257 
258  for(const xAOD::Muon * mu : *muonContainer){
259  const xAOD::TrackParticle* muonTrk = mu->trackParticle(xAOD::Muon::InnerDetectorTrackParticle );
260  if (muonTrk == cascadeVertices[1]->trackParticle(0)) selectedMuons.push_back(mu); //there is always only one muon
261  }
262  ATH_MSG_DEBUG("selectedMuon size "<<selectedMuons.size()); //always only one muon
263 
264  //-----------------------------------
265  // Create link for Muon
266  // create tmp vector of preceding vertex links
267  MuonsLinkVector preMuLinks;
268 
269  // loop over input precedingMuons
270  auto muItr = selectedMuons.begin();
271  for(; muItr != selectedMuons.end(); muItr++) {
272  // create element link
273  MuonsLink muLink;
274  muLink.setElement(*muItr);
276 
277  // sanity check : is the link valid?
278  if( !muLink.isValid() ) continue;
279 
280  // link is OK, store it in the tmp vector
281  preMuLinks.push_back( muLink );
282 
283  } // end of loop over preceding vertices
284 
285  // all OK: store preceding vertex links in the aux store
286  MuonsLinksDecor(*cascadeVertices[1]) = preMuLinks;
287  //-----------------------------------
288  //----------------------------------- mu
289 
290  //----------------------------------- Dx / Lambda_c
291  // Identify the input D_(s)+ or Lambda_c+
292  const xAOD::Vertex* dxVertex = BPhysPVCascadeTools::FindVertex<3>(dxContainer, cascadeVertices[0]);
293  ATH_MSG_DEBUG("1 pt D_(s)+/Lambda_c+ tracks " << cascadeVertices[1]->trackParticle(0)->pt()<<", "
294  << cascadeVertices[0]->trackParticle(0)->pt()<< ", "<< cascadeVertices[0]->trackParticle(1)->pt()<< ", "
295  << cascadeVertices[0]->trackParticle(2)->pt());
296  if (dxVertex) ATH_MSG_DEBUG("2 pt D_(s)+/Lambda_c+ tracks " << dxVertex->trackParticle(0)->pt() << ", "
297  << dxVertex->trackParticle(1)->pt() << ", " << dxVertex->trackParticle(2)->pt());
298  // Set links to input vertices
299  std::vector<const xAOD::Vertex*> dxVerticestoLink;
300  if (dxVertex) dxVerticestoLink.push_back(dxVertex);
301  else ATH_MSG_WARNING("Could not find linking D_(s)+/Lambda_c+");
302  if(!BPhysPVCascadeTools::LinkVertices(DxLinksDecor, dxVerticestoLink, dxContainer, cascadeVertices[1]))
303  ATH_MSG_ERROR("Error decorating with D_(s)+/Lambda_c+ vertices");
304 
305  // tag for D- and Lambda_c-
306  bool tagDp(true);
307  if (dxVertex) {
308  if((std::abs(m_Dx_pid)==411 || std::abs(m_Dx_pid)==4122) && (dxVertex->trackParticle(2)->charge()==-1)) tagDp = false;
309  }
310  //----------------------------------- Dx / Lambda_c
311 
312  double mass_b = m_vtx0MassHypo;
313  double mass_d = m_vtx1MassHypo;
314  std::vector<double> massesMu;
315  massesMu.push_back(m_vtx0Daug1MassHypo); // µ
316  std::vector<double> massesDx;
317  if(tagDp){
318  massesDx.push_back(m_vtx1Daug1MassHypo); // K (Ds+/-) or π (D+, Lambda_c+)
319  massesDx.push_back(m_vtx1Daug2MassHypo); // K
320  }else{ // Change the order for D- or Lambda_c-
321  massesDx.push_back(m_vtx1Daug2MassHypo); // K
322  massesDx.push_back(m_vtx1Daug1MassHypo); // π
323  }
324  massesDx.push_back(m_vtx1Daug3MassHypo); // π (Ds+/-, D+-) or p (Lambda_c+)
325  std::vector<double> Masses;
326  Masses.push_back(m_vtx0Daug1MassHypo); // µ
327  Masses.push_back(m_vtx1MassHypo); // Dx / Lambda_c
328 
329  // loop over candidates -- Don't apply PV_minNTracks requirement here
330  // because it may result in exclusion of the high-pt PV.
331  // get good PVs
332 
333 
334  xAOD::BPhysHypoHelper vtx(m_hypoName, mainVertex);
335 
336  // Get refitted track momenta from all vertices, charged tracks only
338 
339  // Decorate main vertex
340  //
341  // 1.a) mass, mass error
342  BPHYS_CHECK( vtx.setMass(m_CascadeTools->invariantMass(moms[1])) );
343  BPHYS_CHECK( vtx.setMassErr(m_CascadeTools->invariantMassError(moms[1],x->getCovariance()[1])) );
344  // 1.b) pt and pT error (the default pt of mainVertex is != the pt of the full cascade fit!)
345  Pt_decor(*mainVertex) = m_CascadeTools->pT(moms[1]);
346  PtErr_decor(*mainVertex) = m_CascadeTools->pTError(moms[1],x->getCovariance()[1]);
347  // 1.c) chi2 and ndof (the default chi2 of mainVertex is != the chi2 of the full cascade fit!)
348  chi2_decor(*mainVertex) = x->fitChi2();
349  ndof_decor(*mainVertex) = x->nDoF();
350  Dx_chi2_svdecor(*mainVertex) = cascadeVertices[0]->chiSquared();
351  Dx_ndof_svdecor(*mainVertex) = cascadeVertices[0]->numberDoF();
352 
353  //--------------------- mu
354  float muM = 0., mupt = 0., mueta = 0.;
355  if (selectedMuons.size()==1) {
356  TLorentzVector p4_mu1;
357  p4_mu1.SetPtEtaPhiM(selectedMuons.at(0)->pt(),
358  selectedMuons.at(0)->eta(),
359  selectedMuons.at(0)->phi(), m_vtx0Daug1MassHypo);
360  muM = p4_mu1.M();
361  mupt = p4_mu1.Pt();
362  mueta = p4_mu1.Eta();
363  }
364  MuMass_decor(*mainVertex) = muM;
365  MuPt_decor(*mainVertex) = mupt;
366  MuEta_decor(*mainVertex) = mueta;
367  ChargeMu_decor(*mainVertex) = selectedMuons.at(0)->charge();
368  MuChi2_decor(*mainVertex) = cascadeVertices[1]->trackParticle(0)->chiSquared();
369  MunDoF_decor(*mainVertex) = cascadeVertices[1]->trackParticle(0)->numberDoF();
370 
371  // track contribution to chi2 of cascasde fit
372  std::vector< Trk::VxTrackAtVertex > trkAtB = cascadeVertices[1]->vxTrackAtVertex();
373  MuChi2B_decor(*mainVertex) = trkAtB.at(0).trackQuality().chiSquared();
374  MunDoFB_decor(*mainVertex) = trkAtB.at(0).trackQuality().numberDoF();
375 
376  //--------------------- Dx / Lambda_c
377  //tagDp = true by default, for Ds+/-, D+ and Lambda_c+
378  float massKX1 = 0.;
379  float massKX1X2 = 0.;
380  float RapidityKX1X2 = 0.;
381  if (dxVertex) {
382  TLorentzVector p4_h1, p4_h2, p4_h3;
383  if(tagDp){
384  p4_h1.SetPtEtaPhiM(dxVertex->trackParticle(0)->pt(),
385  dxVertex->trackParticle(0)->eta(),
386  dxVertex->trackParticle(0)->phi(), m_vtx1Daug1MassHypo);
387  ChargeX1_decor(*mainVertex) = dxVertex->trackParticle(0)->charge();
388  p4_h2.SetPtEtaPhiM(dxVertex->trackParticle(1)->pt(),
389  dxVertex->trackParticle(1)->eta(),
390  dxVertex->trackParticle(1)->phi(), m_vtx1Daug2MassHypo);
391  ChargeK_decor(*mainVertex) = dxVertex->trackParticle(1)->charge();
392  }else{ // Change the order for D- and Lambda_c-
393  p4_h1.SetPtEtaPhiM(dxVertex->trackParticle(0)->pt(),
394  dxVertex->trackParticle(0)->eta(),
395  dxVertex->trackParticle(0)->phi(), m_vtx1Daug2MassHypo);
396  ChargeK_decor(*mainVertex) = dxVertex->trackParticle(0)->charge();
397 
398  p4_h2.SetPtEtaPhiM(dxVertex->trackParticle(1)->pt(),
399  dxVertex->trackParticle(1)->eta(),
400  dxVertex->trackParticle(1)->phi(), m_vtx1Daug1MassHypo);
401  ChargeX1_decor(*mainVertex) = dxVertex->trackParticle(1)->charge(); // K (Ds+) or π (D+)
402  }
403  p4_h3.SetPtEtaPhiM(dxVertex->trackParticle(2)->pt(),
404  dxVertex->trackParticle(2)->eta(),
405  dxVertex->trackParticle(2)->phi(), m_vtx1Daug3MassHypo);
406  ChargeX2_decor(*mainVertex) = dxVertex->trackParticle(2)->charge(); // π (Ds+, D+) or p (Lambda_c+)
407 
408  massKX1 = (p4_h1 + p4_h2).M();
409  massKX1X2 = (p4_h1 + p4_h2 + p4_h3).M();
410  RapidityKX1X2 = (p4_h1 + p4_h2 + p4_h3).Rapidity();
411  }
412  massKX1_svdecor(*mainVertex) = massKX1;
413  massKX1X2_svdecor(*mainVertex) = massKX1X2;
414  RapidityKX1X2_svdecor(*mainVertex) = RapidityKX1X2;
415 
416 
417  ATH_CHECK(helper.FillCandwithRefittedVertices(m_refitPV, pvContainer,
418  refPvContainer, &(*m_pvRefitter), m_PV_max, m_DoVertexType, x, 1, mass_b, vtx));
419 
420 
421  // 4) decorate the main vertex with V0 vertex mass, pt, lifetime and lxy values (plus errors)
422  // V0 points to the main vertex, so lifetime and lxy are w.r.t the main vertex
423  Mass_svdecor(*mainVertex) = m_CascadeTools->invariantMass(moms[0]);
424  MassErr_svdecor(*mainVertex) = m_CascadeTools->invariantMassError(moms[0],x->getCovariance()[0]);
425  Pt_svdecor(*mainVertex) = m_CascadeTools->pT(moms[0]);
426  PtErr_svdecor(*mainVertex) = m_CascadeTools->pTError(moms[0],x->getCovariance()[0]);
427  Lxy_svdecor(*mainVertex) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
428  LxyErr_svdecor(*mainVertex) = m_CascadeTools->lxyError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
429  Tau_svdecor(*mainVertex) = m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1]);
430  TauErr_svdecor(*mainVertex) = m_CascadeTools->tauError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
431 
432  // Some checks in DEBUG mode
433  ATH_MSG_DEBUG("chi2 " << x->fitChi2()
434  << " chi2_1 " << m_V0Tools->chisq(cascadeVertices[0])
435  << " chi2_2 " << m_V0Tools->chisq(cascadeVertices[1])
436  << " vprob " << m_CascadeTools->vertexProbability(x->nDoF(),x->fitChi2()));
437  ATH_MSG_DEBUG("ndf " << x->nDoF() << " ndf_1 " << m_V0Tools->ndof(cascadeVertices[0]) << " ndf_2 " << m_V0Tools->ndof(cascadeVertices[1]));
438  ATH_MSG_DEBUG("V0Tools mass_d " << m_V0Tools->invariantMass(cascadeVertices[0],massesDx)
439  << " error " << m_V0Tools->invariantMassError(cascadeVertices[0],massesDx)
440  << " mass_J " << m_V0Tools->invariantMass(cascadeVertices[1],massesMu)
441  << " error " << m_V0Tools->invariantMassError(cascadeVertices[1],massesMu));
442  // masses and errors, using track masses assigned in the fit
443  double Mass_B = m_CascadeTools->invariantMass(moms[1]);
444  double Mass_D = m_CascadeTools->invariantMass(moms[0]);
445  double Mass_B_err = m_CascadeTools->invariantMassError(moms[1],x->getCovariance()[1]);
446  double Mass_D_err = m_CascadeTools->invariantMassError(moms[0],x->getCovariance()[0]);
447  ATH_MSG_DEBUG("Mass_B " << Mass_B << " Mass_D " << Mass_D);
448  ATH_MSG_DEBUG("Mass_B_err " << Mass_B_err << " Mass_D_err " << Mass_D_err);
449  double mprob_B = m_CascadeTools->massProbability(mass_b,Mass_B,Mass_B_err);
450  double mprob_D = m_CascadeTools->massProbability(mass_d,Mass_D,Mass_D_err);
451  ATH_MSG_DEBUG("mprob_B " << mprob_B << " mprob_D " << mprob_D);
452  // masses and errors, assigning user defined track masses
453  ATH_MSG_DEBUG("Mass_b " << m_CascadeTools->invariantMass(moms[1],Masses)
454  << " Mass_d " << m_CascadeTools->invariantMass(moms[0],massesDx));
455  ATH_MSG_DEBUG("Mass_b_err " << m_CascadeTools->invariantMassError(moms[1],x->getCovariance()[1],Masses)
456  << " Mass_d_err " << m_CascadeTools->invariantMassError(moms[0],x->getCovariance()[0],massesDx));
457  ATH_MSG_DEBUG("pt_b " << m_CascadeTools->pT(moms[1])
458  << " pt_d " << m_CascadeTools->pT(moms[0])
459  << " pt_dp " << m_V0Tools->pT(cascadeVertices[0]));
460  ATH_MSG_DEBUG("ptErr_b " << m_CascadeTools->pTError(moms[1],x->getCovariance()[1])
461  << " ptErr_d " << m_CascadeTools->pTError(moms[0],x->getCovariance()[0])
462  << " ptErr_dp " << m_V0Tools->pTError(cascadeVertices[0]));
463  ATH_MSG_DEBUG("lxy_B " << m_V0Tools->lxy(cascadeVertices[1],primaryVertex) << " lxy_D " << m_V0Tools->lxy(cascadeVertices[0],cascadeVertices[1]));
464  ATH_MSG_DEBUG("lxy_b " << m_CascadeTools->lxy(moms[1],cascadeVertices[1],primaryVertex) << " lxy_d " << m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]));
465  ATH_MSG_DEBUG("lxyErr_b " << m_CascadeTools->lxyError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
466  << " lxyErr_d " << m_CascadeTools->lxyError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
467  << " lxyErr_dp " << m_V0Tools->lxyError(cascadeVertices[0],cascadeVertices[1]));
468  ATH_MSG_DEBUG("tau_B " << m_CascadeTools->tau(moms[1],cascadeVertices[1],primaryVertex,mass_b)
469  << " tau_dp " << m_V0Tools->tau(cascadeVertices[0],cascadeVertices[1],massesDx));
470  ATH_MSG_DEBUG("tau_b " << m_CascadeTools->tau(moms[1],cascadeVertices[1],primaryVertex)
471  << " tau_d " << m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1])
472  << " tau_D " << m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1],mass_d));
473  ATH_MSG_DEBUG("tauErr_b " << m_CascadeTools->tauError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
474  << " tauErr_d " << m_CascadeTools->tauError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
475  << " tauErr_dp " << m_V0Tools->tauError(cascadeVertices[0],cascadeVertices[1],massesDx));
476  ATH_MSG_DEBUG("TauErr_b " << m_CascadeTools->tauError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex,mass_b)
477  << " TauErr_d " << m_CascadeTools->tauError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1],mass_d)
478  << " TauErr_dp " << m_V0Tools->tauError(cascadeVertices[0],cascadeVertices[1],massesDx,mass_d));
479 
480  ATH_MSG_DEBUG("CascadeTools main vert wrt PV " << " CascadeTools SV " << " V0Tools SV");
481  ATH_MSG_DEBUG("a0z " << m_CascadeTools->a0z(moms[1],cascadeVertices[1],primaryVertex)
482  << ", " << m_CascadeTools->a0z(moms[0],cascadeVertices[0],cascadeVertices[1])
483  << ", " << m_V0Tools->a0z(cascadeVertices[0],cascadeVertices[1]));
484  ATH_MSG_DEBUG("a0zErr " << m_CascadeTools->a0zError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
485  << ", " << m_CascadeTools->a0zError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
486  << ", " << m_V0Tools->a0zError(cascadeVertices[0],cascadeVertices[1]));
487  ATH_MSG_DEBUG("a0xy " << m_CascadeTools->a0xy(moms[1],cascadeVertices[1],primaryVertex)
488  << ", " << m_CascadeTools->a0xy(moms[0],cascadeVertices[0],cascadeVertices[1])
489  << ", " << m_V0Tools->a0xy(cascadeVertices[0],cascadeVertices[1]));
490  ATH_MSG_DEBUG("a0xyErr " << m_CascadeTools->a0xyError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
491  << ", " << m_CascadeTools->a0xyError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
492  << ", " << m_V0Tools->a0xyError(cascadeVertices[0],cascadeVertices[1]));
493  ATH_MSG_DEBUG("a0 " << m_CascadeTools->a0(moms[1],cascadeVertices[1],primaryVertex)
494  << ", " << m_CascadeTools->a0(moms[0],cascadeVertices[0],cascadeVertices[1])
495  << ", " << m_V0Tools->a0(cascadeVertices[0],cascadeVertices[1]));
496  ATH_MSG_DEBUG("a0Err " << m_CascadeTools->a0Error(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
497  << ", " << m_CascadeTools->a0Error(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
498  << ", " << m_V0Tools->a0Error(cascadeVertices[0],cascadeVertices[1]));
499  ATH_MSG_DEBUG("x0 " << m_V0Tools->vtx(cascadeVertices[0]).x() << " y0 " << m_V0Tools->vtx(cascadeVertices[0]).y() << " z0 " << m_V0Tools->vtx(cascadeVertices[0]).z());
500  ATH_MSG_DEBUG("x1 " << m_V0Tools->vtx(cascadeVertices[1]).x() << " y1 " << m_V0Tools->vtx(cascadeVertices[1]).y() << " z1 " << m_V0Tools->vtx(cascadeVertices[1]).z());
501  ATH_MSG_DEBUG("X0 " << primaryVertex->x() << " Y0 " << primaryVertex->y() << " Z0 " << primaryVertex->z());
502  ATH_MSG_DEBUG("rxy0 " << m_V0Tools->rxy(cascadeVertices[0]) << " rxyErr0 " << m_V0Tools->rxyError(cascadeVertices[0]));
503  ATH_MSG_DEBUG("rxy1 " << m_V0Tools->rxy(cascadeVertices[1]) << " rxyErr1 " << m_V0Tools->rxyError(cascadeVertices[1]));
504  ATH_MSG_DEBUG("Rxy0 wrt PV " << m_V0Tools->rxy(cascadeVertices[0],primaryVertex) << " RxyErr0 wrt PV " << m_V0Tools->rxyError(cascadeVertices[0],primaryVertex));
505  ATH_MSG_DEBUG("Rxy1 wrt PV " << m_V0Tools->rxy(cascadeVertices[1],primaryVertex) << " RxyErr1 wrt PV " << m_V0Tools->rxyError(cascadeVertices[1],primaryVertex));
506  ATH_MSG_DEBUG("number of covariance matrices " << (x->getCovariance()).size());
507 
508  } // loop over cascadeinfoContainer
509 
510  // Deleting cascadeinfo since this won't be stored.
511  // Vertices have been kept in m_cascadeOutputs and should be owned by their container
512  for (auto x : cascadeinfoContainer) delete x;
513 
514  return StatusCode::SUCCESS;
515  }
516 
517 
518  MuPlusDsCascade::MuPlusDsCascade(const std::string& t, const std::string& n, const IInterface* p) : AthAlgTool(t,n,p),
519  m_vertexDxContainerKey(""),
520  m_cascadeOutputsKeys{ "MuPlusDsCascadeVtx1", "MuPlusDsCascadeVtx2" },
521  m_VxPrimaryCandidateName("PrimaryVertices"),
522  m_DxMassLower(0.0),
523  m_DxMassUpper(10000.0),
524  m_MassLower(0.0),
525  m_MassUpper(20000.0),
526  m_vtx0MassHypo(-1),
527  m_vtx1MassHypo(-1),
528  m_vtx0Daug1MassHypo(-1),
529  m_vtx1Daug1MassHypo(-1),
530  m_vtx1Daug2MassHypo(-1),
531  m_vtx1Daug3MassHypo(-1),
532  m_Dx_pid(431),
533  m_constrDx(true),
534  m_chi2cut(-1.0),
535  m_iVertexFitter("Trk::TrkVKalVrtFitter"),
536  m_pvRefitter("Analysis::PrimaryVertexRefitter",this),
537  m_V0Tools("Trk::V0Tools"),
538  m_CascadeTools("DerivationFramework::CascadeTools"),
539  m_muonCollectionKey("StacoMuonCollection"),
540  m_trkSelector("InDet::TrackSelectorTool"),
541  m_thresholdPt(0.0),
542  m_mcpCuts(true),
543  m_combOnly(false),
544  m_useCombMeasurement(false)
545  {
546  declareProperty("DxVertices", m_vertexDxContainerKey);
547  declareProperty("VxPrimaryCandidateName", m_VxPrimaryCandidateName);
548  declareProperty("RefPVContainerName", m_refPVContainerName = "RefittedPrimaryVertices");
549  declareProperty("DxMassLowerCut", m_DxMassLower);
550  declareProperty("DxMassUpperCut", m_DxMassUpper);
551  declareProperty("MassLowerCut", m_MassLower);
552  declareProperty("MassUpperCut", m_MassUpper);
553  declareProperty("HypothesisName", m_hypoName = "B");
554  declareProperty("Vtx0MassHypo", m_vtx0MassHypo);
555  declareProperty("Vtx1MassHypo", m_vtx1MassHypo);
556  declareProperty("Vtx0Daug1MassHypo", m_vtx0Daug1MassHypo);
557  declareProperty("Vtx1Daug1MassHypo", m_vtx1Daug1MassHypo);
558  declareProperty("Vtx1Daug2MassHypo", m_vtx1Daug2MassHypo);
559  declareProperty("Vtx1Daug3MassHypo", m_vtx1Daug3MassHypo);
560  declareProperty("DxHypothesis", m_Dx_pid);
561  declareProperty("ApplyDxMassConstraint", m_constrDx);
562  declareProperty("Chi2Cut", m_chi2cut);
563  declareProperty("RefitPV", m_refitPV = true);
564  declareProperty("MaxnPV", m_PV_max = 999);
565  declareProperty("MinNTracksInPV", m_PV_minNTracks = 0);
566  declareProperty("DoVertexType", m_DoVertexType = 7);
567  declareProperty("TrkVertexFitterTool", m_iVertexFitter);
568  declareProperty("PVRefitter", m_pvRefitter);
569  declareProperty("V0Tools", m_V0Tools);
570  declareProperty("CascadeTools", m_CascadeTools);
571  declareProperty("CascadeVertexCollections", m_cascadeOutputsKeys);
572  declareProperty("muonCollectionKey", m_muonCollectionKey);
573  declareProperty("TrackSelectorTool", m_trkSelector);
574  declareProperty("muonThresholdPt", m_thresholdPt);
575  declareProperty("useMCPCuts", m_mcpCuts);
576  declareProperty("combOnly", m_combOnly);
577  declareProperty("useCombinedMeasurement", m_useCombMeasurement);
578  }
579 
581 
582  StatusCode MuPlusDsCascade::performSearch(std::vector<Trk::VxCascadeInfo*> *cascadeinfoContainer) const
583  {
584  ATH_MSG_DEBUG( "MuPlusDsCascade::performSearch" );
585  assert(cascadeinfoContainer!=nullptr);
586 
587  // Get TrackParticle container (for setting links to the original tracks)
588  const xAOD::TrackParticleContainer *trackContainer(nullptr);
589  ATH_CHECK(evtStore()->retrieve(trackContainer , "InDetTrackParticles" ));
590 
591  // Get V0 container
592  const xAOD::VertexContainer *dxContainer(nullptr);
594 
595  double mass_d = m_vtx1MassHypo;
596  std::vector<const xAOD::TrackParticle*> tracksMu;
597  std::vector<const xAOD::TrackParticle*> tracksDx;
598  std::vector<double> massesMu;
599  massesMu.push_back(m_vtx0Daug1MassHypo); // µ
600  std::vector<double> massesDx;
601  massesDx.push_back(m_vtx1Daug1MassHypo); // K (Ds+/-) or π (D+, Lambda_c+)
602  massesDx.push_back(m_vtx1Daug2MassHypo); // K
603  massesDx.push_back(m_vtx1Daug3MassHypo); // π (Ds+/-, D+) ot p (Lambda_c+)
604  std::vector<double> massesDm; // Alter the order of masses for D- and Lambda_c-
605  massesDm.push_back(m_vtx1Daug2MassHypo); // K
606  massesDm.push_back(m_vtx1Daug1MassHypo); // π
607  massesDm.push_back(m_vtx1Daug3MassHypo); // π (D-) or p (Lambda_c-)
608  std::vector<double> Masses;
609  Masses.push_back(m_vtx0Daug1MassHypo); // µ
610  Masses.push_back(m_vtx1MassHypo); // Dx / Lambda_c
611 
612  //-------------------------------------------------------------------------------
613  // Retrieving and selecting muons (taken from JpsiFinder.cxx)
614  // Get the muons from StoreGate
615  const xAOD::MuonContainer* importedMuonCollection;
616  StatusCode sc = evtStore()->retrieve(importedMuonCollection,m_muonCollectionKey);
617 
618  if(sc.isFailure()){
619  ATH_MSG_WARNING("No muon collection with key " << m_muonCollectionKey << " found in StoreGate");
620  return StatusCode::SUCCESS;;
621  }else{
622  ATH_MSG_DEBUG("Found muon collections with key "<<m_muonCollectionKey);
623  }
624  ATH_MSG_DEBUG("Muon container size "<<importedMuonCollection->size());
625 
626  // Typedef for vectors of muons
627  typedef std::vector<const xAOD::Muon*> MuonBag;
628 
629  // Select the muons
630  const xAOD::Vertex* vx = 0;
631  MuonBag theMuonsAfterSelection;
632  for (auto mu : *importedMuonCollection) {
633  if ( mu == nullptr ) continue;
634  if (!mu->inDetTrackParticleLink().isValid()) continue; // No muons without ID tracks
635  const xAOD::TrackParticle* muonTrk = *(mu->inDetTrackParticleLink());
636  if ( muonTrk==nullptr) continue;
637  if ( !m_trkSelector->decision(*muonTrk, vx) ) continue; // all ID tracks must pass basic tracking cuts
638  if ( std::fabs(muonTrk->pt())<m_thresholdPt ) continue; // higher pt cut if needed
639  if ( m_mcpCuts && !mu->passesIDCuts()) continue; // cuts of the MCP group recommendation
640  if ( m_combOnly && mu->muonType() != xAOD::Muon::Combined ) continue; // require combined muons
641  if ( mu->muonType() == xAOD::Muon::SiliconAssociatedForwardMuon && !m_useCombMeasurement) continue;
642 
643  theMuonsAfterSelection.push_back(mu);
644  }
645  if (theMuonsAfterSelection.size() == 0) return StatusCode::SUCCESS;;
646  ATH_MSG_DEBUG("Number of muons after selection: " << theMuonsAfterSelection.size());
647  //-------------------------------------------------------------------------------
648 
649  //-------------------------------------------------------------------------------
650  // Select the D_s+/D+/Lambda_c+ candidates before calling cascade fit
651  std::vector<const xAOD::Vertex*> selectedDxCandidates;
652 
653 
654  for(auto vxcItr : *dxContainer){
655 
656  // Check the passed flag first
657  const xAOD::Vertex* vtx = vxcItr;
658  if(std::abs(m_Dx_pid)==431) { // D_s+/-
659  SG::AuxElement::Accessor<Char_t> flagAcc1("passed_Ds");
660  if(flagAcc1.isAvailable(*vtx)){
661  if(!flagAcc1(*vtx)) continue;
662  }
663  }
664 
665  if(std::abs(m_Dx_pid)==411) { // D+/-
666  SG::AuxElement::Accessor<Char_t> flagAcc1("passed_Dp");
667  SG::AuxElement::Accessor<Char_t> flagAcc2("passed_Dm");
668  bool isDp(true);
669  bool isDm(true);
670  if(flagAcc1.isAvailable(*vtx)){
671  if(!flagAcc1(*vtx)) isDp = false;
672  }
673  if(flagAcc2.isAvailable(*vtx)){
674  if(!flagAcc2(*vtx)) isDm = false;
675  }
676  if(!(isDp||isDm)) continue;
677  }
678 
679  // Track selection - Loose
680  if ( !m_trackSelectionTools->accept(vxcItr->trackParticle(0)) ){
681  ATH_MSG_DEBUG(" Original Dx/Lambda_c candidate rejected by the track's cut level - loose ");
682  continue;
683  }
684  if ( !m_trackSelectionTools->accept(vxcItr->trackParticle(1)) ){
685  ATH_MSG_DEBUG(" Original Dx/Lambda_c candidate rejected by the track's cut level - loose ");
686  continue;
687  }
688  if ( !m_trackSelectionTools->accept(vxcItr->trackParticle(2)) ){
689  ATH_MSG_DEBUG(" Original Dx/Lambda_c candidate rejected by the track's cut level - loose ");
690  continue;
691  }
692 
693  // Ensure the total charge is correct
694  if( std::abs( vxcItr->trackParticle(0)->charge()+vxcItr->trackParticle(1)->charge()+vxcItr->trackParticle(2)->charge() ) != 1 ){
695  ATH_MSG_DEBUG(" Original Dx/Lambda_c candidate rejected by the charge requirement: "
696  << vxcItr->trackParticle(0)->charge() << ", " << vxcItr->trackParticle(1)->charge() << ", " << vxcItr->trackParticle(2)->charge() );
697  continue;
698  }
699 
700  // Check D_(s)/Lambda_c +/- candidate invariant mass and skip if need be
701  double mass_D;
702  if( (std::abs(m_Dx_pid) == 411 || std::abs(m_Dx_pid) == 4122) && vxcItr->trackParticle(2)->charge()<0) // D- & Lambda_c-
703  mass_D = m_V0Tools->invariantMass(vxcItr,massesDm);
704  else // D+, D_s+/-, Lambda_c+
705  mass_D = m_V0Tools->invariantMass(vxcItr,massesDx);
706  ATH_MSG_DEBUG("D_(s)/Lambda_c mass " << mass_D);
707  if(mass_D < m_DxMassLower || mass_D > m_DxMassUpper) {
708  ATH_MSG_DEBUG(" Original D_(s)/Lambda_c candidate rejected by the mass cut: mass = "
709  << mass_D << " != (" << m_DxMassLower << ", " << m_DxMassUpper << ")" );
710  continue;
711  }
712 
713  // Add loose cut on phi(K+K-) mass for D_s->phi π
714  if(std::abs(m_Dx_pid)==431){
715  TLorentzVector p4Kp_in, p4Km_in;
716  p4Kp_in.SetPtEtaPhiM( vxcItr->trackParticle(0)->pt(),
717  vxcItr->trackParticle(0)->eta(),
718  vxcItr->trackParticle(0)->phi(), m_vtx1Daug1MassHypo);
719  p4Km_in.SetPtEtaPhiM( vxcItr->trackParticle(1)->pt(),
720  vxcItr->trackParticle(1)->eta(),
721  vxcItr->trackParticle(1)->phi(), m_vtx1Daug2MassHypo);
722  double mass_phi = (p4Kp_in + p4Km_in).M();
723  ATH_MSG_DEBUG("phi mass " << mass_phi);
724  if(mass_phi > 1600) { //1st cut on phi(K+K-) mass, loose one
725  ATH_MSG_DEBUG(" Original phi candidate rejected by the mass cut: mass = " << mass_phi );
726  continue;
727 
728  }
729  }
730 
731  // Cut on pT(p) for Lambda_c
732  if(std::abs(m_Dx_pid)==4122){
733  if (vxcItr->trackParticle(2)->pt() < 2400) continue;
734  }
735 
736  // L.G., 26.03.21: Check there are no previous D+/- combinations with inverse pion positions
737  if(std::abs(m_Dx_pid)==411 && selectedDxCandidates.size()>0) {
738  bool Dpmcopy(false);
739  for(auto dxItr : selectedDxCandidates){
740  if(vxcItr->trackParticle(2)/*π*/->charge()>0) { // D+
741  if ( vxcItr->trackParticle(2) == dxItr->trackParticle(0) && vxcItr->trackParticle(0) == dxItr->trackParticle(2)
742  && vxcItr->trackParticle(1) == dxItr->trackParticle(1) ) Dpmcopy = true;
743  } else { // D-
744  if ( vxcItr->trackParticle(2) == dxItr->trackParticle(1) && vxcItr->trackParticle(1) == dxItr->trackParticle(2)
745  && vxcItr->trackParticle(0) == dxItr->trackParticle(0) ) Dpmcopy = true;
746  }
747  }
748  if (Dpmcopy) continue;
749  }
750 
751  selectedDxCandidates.push_back(vxcItr);
752  } // end for(auto vxcItr : dxContainer)
753  if(selectedDxCandidates.size()<1) return StatusCode::SUCCESS;
754  //-------------------------------------------------------------------------------
755 
756  //-------------------------------------------------------------------------------
757  // Select mu D_(s)+/Lambda_c candidates
758  // Iterate over muons
759  for(auto muItr : theMuonsAfterSelection){
760  tracksMu.clear();
761  //Convert to trackParticle base
762  auto TrkMuon = muItr->trackParticle( xAOD::Muon::InnerDetectorTrackParticle );
763  tracksMu.push_back(TrkMuon);
764  if (tracksMu.size() != 1 || massesMu.size() != 1 ) {
765  ATH_MSG_WARNING("Problems with muon input");
766  }
767 
768  // Iterate over D_(s)/Lambda_c +/- vertices
769  for(auto dxItr : selectedDxCandidates){
770  // Check identical tracks in input
771  if(std::find(tracksMu.cbegin(), tracksMu.cend(), dxItr->trackParticle(0)) != tracksMu.cend()) continue;
772  if(std::find(tracksMu.cbegin(), tracksMu.cend(), dxItr->trackParticle(1)) != tracksMu.cend()) continue;
773  if(std::find(tracksMu.cbegin(), tracksMu.cend(), dxItr->trackParticle(2)) != tracksMu.cend()) continue;
774 
775  size_t dxTrkNum = dxItr->nTrackParticles();
776  tracksDx.clear();
777  for( unsigned int it=0; it<dxTrkNum; it++) tracksDx.push_back(dxItr->trackParticle(it));
778  if (tracksDx.size() != 3 || massesDx.size() != 3 ) {
779  ATH_MSG_WARNING("Problems with D_(s)/Lambda_c +/- input");
780  }
781 
782  if(std::find(tracksDx.cbegin(), tracksDx.cend(), TrkMuon) != tracksDx.cend()) continue; //looking for the muons in the D+ tracks
783 
784  ATH_MSG_DEBUG("Using tracks" << tracksMu[0] << ", " << tracksDx[0] << ", " << tracksDx[1] << ", " << tracksDx[2]);
785  // Apply the user's settings to the fitter
786  // Reset
787  std::unique_ptr<Trk::IVKalState> state (m_iVertexFitter->makeState());
788  // Robustness
789  int robustness = 0;
790  m_iVertexFitter->setRobustness(robustness, *state);
791  // Build up the topology
792  // Vertex list
793  std::vector<Trk::VertexID> vrtList;
794  // D_(s)/Lambda_c +/- vertex
795  Trk::VertexID vID;
796  if (m_constrDx) {
797  if((std::abs(m_Dx_pid)==411 || std::abs(m_Dx_pid)==4122 ) && dxItr->trackParticle(2)->charge()<0) // D-, Lambda_c-
798  vID = m_iVertexFitter->startVertex(tracksDx,massesDm, *state,mass_d);
799  else // D+, D_s+/-, Lambda_c+
800  vID = m_iVertexFitter->startVertex(tracksDx,massesDx, *state,mass_d);
801  } else {
802  if((std::abs(m_Dx_pid)==411 || std::abs(m_Dx_pid)==4122 ) && dxItr->trackParticle(2)->charge()<0) // D-, Lambda_c-
803  vID = m_iVertexFitter->startVertex(tracksDx,massesDm, *state);
804  else // D+, D_s+/-, Lambda_c+
805  vID = m_iVertexFitter->startVertex(tracksDx,massesDx, *state);
806  }
807  vrtList.push_back(vID);
808 
809  // B vertex including muon
810  //For mass constraint, not used yet
811  Trk::VertexID vID2 = m_iVertexFitter->nextVertex(tracksMu,massesMu,vrtList, *state);
812  ATH_MSG_DEBUG(vID2);
813 
814  // Do the work
815  std::unique_ptr<Trk::VxCascadeInfo> result(m_iVertexFitter->fitCascade(*state));
816 
817  if (result != nullptr) {
818  // reset links to original tracks
819  BPhysPVCascadeTools::PrepareVertexLinks(result.get(), trackContainer);
820  ATH_MSG_DEBUG("storing tracks " << ((result->vertices())[0] /*D_(s)+*/)->trackParticle(0) << ", "
821  << ((result->vertices())[0])->trackParticle(1) << ", "
822  << ((result->vertices())[0])->trackParticle(2) << ", "
823  << ((result->vertices())[1]/*mu*/)->trackParticle(0));
824 
825  // necessary to prevent memory leak
826  result->setSVOwnership(true);
827 
828  // Chi2/DOF cut
829  double bChi2DOF = result->fitChi2()/result->nDoF();
830  ATH_MSG_DEBUG("Candidate chi2/DOF is " << bChi2DOF);
831  bool chi2CutPassed = (m_chi2cut <= 0.0 || bChi2DOF < m_chi2cut);
832 
833  const std::vector<xAOD::Vertex*> &cascadeVertices = result->vertices();
834  const std::vector< std::vector<TLorentzVector> > &moms = result->getParticleMoms();
835 
836  //----------------------------------------------------
837  // retrieve primary vertices
838  //----------------------------------------------------
839  const xAOD::Vertex * primaryVertex(nullptr);
840  const xAOD::VertexContainer *pvContainer(nullptr);
842  ATH_MSG_DEBUG("Found " << m_VxPrimaryCandidateName << " in StoreGate!");
843 
844  if (pvContainer->size()==0){
845  ATH_MSG_WARNING("You have no primary vertices: " << pvContainer->size());
846  return StatusCode::RECOVERABLE;
847  } else {
848  primaryVertex = (*pvContainer)[0];
849  }
850 
851  // Add stronger cut on phi(K+K-) mass for D_s->phi π after cascade fit
852  if(std::abs(m_Dx_pid)==431){
853  TLorentzVector p4Kp_in, p4Km_in;
854  p4Kp_in.SetPtEtaPhiM( cascadeVertices[0]->trackParticle(0)->pt(),
855  cascadeVertices[0]->trackParticle(0)->eta(),
856  cascadeVertices[0]->trackParticle(0)->phi(), m_vtx1Daug1MassHypo);
857  p4Km_in.SetPtEtaPhiM( cascadeVertices[0]->trackParticle(1)->pt(),
858  cascadeVertices[0]->trackParticle(1)->eta(),
859  cascadeVertices[0]->trackParticle(1)->phi(), m_vtx1Daug2MassHypo);
860  double mass_phi = (p4Kp_in + p4Km_in).M();
861  ATH_MSG_DEBUG("phi mass " << mass_phi);
862  if(mass_phi > 1100) {
863  ATH_MSG_DEBUG(" Original phi candidate rejected by the mass cut: mass = " << mass_phi );
864  continue;
865  }
866  }
867 
868  double mass = m_CascadeTools->invariantMass(moms[1]);
869  if(chi2CutPassed) {
870  if (mass >= m_MassLower && mass <= m_MassUpper) {
871  if (m_CascadeTools->lxy(moms[1],cascadeVertices[1],primaryVertex) > 0.09){ //B_Lxy
872  if (m_CascadeTools->pT(moms[1]) > 9500){ //B_pT
873  if (m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]) > 0){ //Dx_lxy
874  if (std::abs(m_Dx_pid)!=411) cascadeinfoContainer->push_back(result.release());
875  else if (m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]) > 0.09){ //D+_Lxy
876  cascadeinfoContainer->push_back(result.release());
877  }//D+_Lxy
878  }//Dx_lxy
879  } //B_pT
880  }//B_Lxy
881  } else {
882  ATH_MSG_DEBUG("Candidate rejected by the mass cut: mass = "
883  << mass << " != (" << m_MassLower << ", " << m_MassUpper << ")" );
884  }
885  } // chi2CutPassed
886 
887  } // result != nullptr
888 
889  } //Iterate over D_(s)/Lambda_c +/- vertices
890 
891  } //Iterate over muons
892 
893  ATH_MSG_DEBUG("cascadeinfoContainer size " << cascadeinfoContainer->size());
894 
895  return StatusCode::SUCCESS;
896  }
897 
898 }
899 
900 
muonContainer
xAOD::MuonContainer * muonContainer
Definition: TrigGlobEffCorrValidation.cxx:188
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
xAOD::BPhysHypoHelper::setMass
bool setMass(const float val)
Set given invariant mass and its error.
Definition: BPhysHypoHelper.cxx:49
xAOD::TrackParticle_v1::pt
virtual double pt() const override final
The transverse momentum ( ) of the particle.
Definition: TrackParticle_v1.cxx:73
DerivationFramework::MuPlusDsCascade::m_partPropSvc
ServiceHandle< IPartPropSvc > m_partPropSvc
Definition: MuPlusDsCascade.h:83
xAOD::Vertex_v1::x
float x() const
Returns the x position.
V0Tools.h
DerivationFramework::MuPlusDsCascade::m_VxPrimaryCandidateName
std::string m_VxPrimaryCandidateName
Name of primary vertex container.
Definition: MuPlusDsCascade.h:60
DerivationFramework::MuPlusDsCascade::MuPlusDsCascade
MuPlusDsCascade(const std::string &t, const std::string &n, const IInterface *p)
Definition: MuPlusDsCascade.cxx:518
DerivationFramework::MuonsLink
ElementLink< xAOD::MuonContainer > MuonsLink
Definition: MuPlusDsCascade.cxx:35
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
DerivationFramework::MuPlusDsCascade::m_eventInfo_key
SG::ReadHandleKey< xAOD::EventInfo > m_eventInfo_key
Definition: MuPlusDsCascade.h:77
DerivationFramework::MuPlusDsCascade::m_vertexDxContainerKey
std::string m_vertexDxContainerKey
Definition: MuPlusDsCascade.h:57
get_generator_info.result
result
Definition: get_generator_info.py:21
DerivationFramework::MuPlusDsCascade::m_muonCollectionKey
std::string m_muonCollectionKey
Definition: MuPlusDsCascade.h:94
xAOD::VertexAuxContainer_v1
Temporary container used until we have I/O for AuxStoreInternal.
Definition: VertexAuxContainer_v1.h:32
Trk::VxCascadeInfo
Definition: VxCascadeInfo.h:75
Trk::VertexID
int VertexID
Definition: IVertexCascadeFitter.h:23
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
DerivationFramework::MuPlusDsCascade::m_combOnly
bool m_combOnly
Definition: MuPlusDsCascade.h:98
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
DerivationFramework::BPhysPVCascadeTools::LinkVertices
static bool LinkVertices(SG::AuxElement::Decorator< VertexLinkVector > &decor, const std::vector< const xAOD::Vertex * > &vertices, const xAOD::VertexContainer *vertexContainer, const xAOD::Vertex *vert)
Definition: BPhysPVCascadeTools.cxx:460
Base_Fragment.mass
mass
Definition: Sherpa_i/share/common/Base_Fragment.py:59
SG::Accessor
Helper class to provide type-safe access to aux data.
Definition: Control/AthContainers/AthContainers/Accessor.h:68
DerivationFramework::MuPlusDsCascade::m_vtx0MassHypo
double m_vtx0MassHypo
Definition: MuPlusDsCascade.h:66
xAOD::TrackParticle_v1::charge
float charge() const
Returns the charge.
Definition: TrackParticle_v1.cxx:150
ConstDataVector.h
DataVector adapter that acts like it holds const pointers.
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:70
xAOD::TrackParticle_v1::eta
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
Definition: TrackParticle_v1.cxx:77
DerivationFramework::BPhysPVCascadeTools::getParticleMass
static double getParticleMass(const HepPDT::ParticleDataTable *pdt, int pdg)
Definition: BPhysPVCascadeTools.cxx:491
DerivationFramework::MuPlusDsCascade::m_PV_minNTracks
size_t m_PV_minNTracks
Definition: MuPlusDsCascade.h:92
DerivationFramework::MuPlusDsCascade::m_CascadeTools
ToolHandle< DerivationFramework::CascadeTools > m_CascadeTools
Definition: MuPlusDsCascade.h:81
DerivationFramework::MuPlusDsCascade::m_V0Tools
ToolHandle< Trk::V0Tools > m_V0Tools
Definition: MuPlusDsCascade.h:80
AuxContainerBase.h
skel.it
it
Definition: skel.GENtoEVGEN.py:396
test_pyathena.pt
pt
Definition: test_pyathena.py:11
ANA_CHECK
#define ANA_CHECK(EXP)
check whether the given expression was successful
Definition: Control/AthToolSupport/AsgMessaging/AsgMessaging/MessageCheck.h:324
LArG4FSStartPointFilter.evt
evt
Definition: LArG4FSStartPointFilter.py:42
DerivationFramework::MuPlusDsCascade::m_trackSelectionTools
std::unique_ptr< InDet::InDetTrackSelectionTool > m_trackSelectionTools
Definition: MuPlusDsCascade.h:82
DerivationFramework::MuPlusDsCascade::m_chi2cut
double m_chi2cut
Definition: MuPlusDsCascade.h:75
DerivationFramework::MuPlusDsCascade::m_MassLower
double m_MassLower
Definition: MuPlusDsCascade.h:64
DerivationFramework::MuPlusDsCascade::m_vtx1Daug3MassHypo
double m_vtx1Daug3MassHypo
Definition: MuPlusDsCascade.h:71
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
xAOD::VertexContainer
VertexContainer_v1 VertexContainer
Definition of the current "Vertex container version".
Definition: VertexContainer.h:14
x
#define x
DerivationFramework::MuonBag
std::vector< const xAOD::Muon * > MuonBag
Definition: BPhysAddMuonBasedInvMass.h:33
xAOD::Muon_v1
Class describing a Muon.
Definition: Muon_v1.h:38
TrkVKalVrtFitter.h
AthenaPoolTestRead.sc
sc
Definition: AthenaPoolTestRead.py:27
runBeamSpotCalibration.helper
helper
Definition: runBeamSpotCalibration.py:112
DerivationFramework::MuPlusDsCascade::m_constrDx
bool m_constrDx
Definition: MuPlusDsCascade.h:74
AthCommonDataStore< AthCommonMsg< AlgTool > >::evtStore
ServiceHandle< StoreGateSvc > & evtStore()
The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.
Definition: AthCommonDataStore.h:85
xAOD::BPhysHypoHelper
Definition: BPhysHypoHelper.h:73
DerivationFramework::BPhysPVCascadeTools::SetVectorInfo
static void SetVectorInfo(xAOD::BPhysHelper &, const Trk::VxCascadeInfo *)
Definition: BPhysPVCascadeTools.cxx:424
BPHYS_CHECK
#define BPHYS_CHECK(EXP)
Useful CHECK macro.
Definition: BPhysHelper.h:738
DerivationFramework::MuPlusDsCascade::initialize
virtual StatusCode initialize() override
Definition: MuPlusDsCascade.cxx:38
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
SG::Decorator
Helper class to provide type-safe access to aux data.
Definition: Decorator.h:59
xAOD::VertexAuxContainer
VertexAuxContainer_v1 VertexAuxContainer
Definition of the current jet auxiliary container.
Definition: VertexAuxContainer.h:19
DerivationFramework::MuPlusDsCascade::m_PV_max
int m_PV_max
Definition: MuPlusDsCascade.h:90
lumiFormat.i
int i
Definition: lumiFormat.py:85
beamspotman.n
n
Definition: beamspotman.py:731
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
BPhysHypoHelper.h
: B-physcis xAOD helpers.
DerivationFramework::MuPlusDsCascade::~MuPlusDsCascade
~MuPlusDsCascade()
Definition: MuPlusDsCascade.cxx:580
DerivationFramework::MuPlusDsCascade::m_DoVertexType
int m_DoVertexType
Definition: MuPlusDsCascade.h:91
DerivationFramework::MuPlusDsCascade::m_pvRefitter
ToolHandle< Analysis::PrimaryVertexRefitter > m_pvRefitter
Definition: MuPlusDsCascade.h:79
DerivationFramework::MuPlusDsCascade::m_mcpCuts
bool m_mcpCuts
Definition: MuPlusDsCascade.h:97
DerivationFramework::MuonsLinkVector
std::vector< MuonsLink > MuonsLinkVector
Definition: MuPlusDsCascade.cxx:36
InDetTrackSelectionTool.h
DerivationFramework::BPhysPVCascadeTools::PrepareVertexLinks
static void PrepareVertexLinks(Trk::VxCascadeInfo *result, const xAOD::TrackParticleContainer *importedTrackCollection)
Definition: BPhysPVCascadeTools.cxx:204
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
DerivationFramework::MuPlusDsCascade::m_vtx0Daug1MassHypo
double m_vtx0Daug1MassHypo
Definition: MuPlusDsCascade.h:68
xAOD::Vertex_v1::trackParticle
const TrackParticle * trackParticle(size_t i) const
Get the pointer to a given track that was used in vertex reco.
Definition: Vertex_v1.cxx:249
DerivationFramework::TrackBag
std::vector< const xAOD::TrackParticle * > TrackBag
Definition: BPhysAddMuonBasedInvMass.h:32
xAOD::Vertex_v1::z
float z() const
Returns the z position.
DerivationFramework
THE reconstruction tool.
Definition: ParticleSortingAlg.h:24
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
DerivationFramework::VertexLink
ElementLink< xAOD::VertexContainer > VertexLink
Definition: Cascade3Plus1.cxx:23
DataVector
Derived DataVector<T>.
Definition: DataVector.h:581
DerivationFramework::MuPlusDsCascade::m_useCombMeasurement
bool m_useCombMeasurement
Definition: MuPlusDsCascade.h:99
DerivationFramework::MuPlusDsCascade::m_cascadeOutputsKeys
std::vector< std::string > m_cascadeOutputsKeys
Definition: MuPlusDsCascade.h:58
BPhysPVCascadeTools.h
DerivationFramework::MuPlusDsCascade::m_vtx1Daug2MassHypo
double m_vtx1Daug2MassHypo
Definition: MuPlusDsCascade.h:70
DerivationFramework::MuPlusDsCascade::m_vtx1MassHypo
double m_vtx1MassHypo
Definition: MuPlusDsCascade.h:67
Trk::Combined
@ Combined
Definition: TrackSummaryTool.h:32
DerivationFramework::MuPlusDsCascade::m_thresholdPt
double m_thresholdPt
Definition: MuPlusDsCascade.h:96
MuPlusDsCascade.h
DerivationFramework::MuPlusDsCascade::m_DxMassUpper
double m_DxMassUpper
Definition: MuPlusDsCascade.h:63
DerivationFramework::MuPlusDsCascade::m_refitPV
bool m_refitPV
Definition: MuPlusDsCascade.h:85
DerivationFramework::BPhysPVCascadeTools
Definition: BPhysPVCascadeTools.h:34
CascadeTools.h
DerivationFramework::MuPlusDsCascade::m_trkSelector
ToolHandle< Trk::ITrackSelectorTool > m_trkSelector
Definition: MuPlusDsCascade.h:95
DerivationFramework::MuPlusDsCascade::performSearch
StatusCode performSearch(std::vector< Trk::VxCascadeInfo * > *cascadeinfoContainer) const
Definition: MuPlusDsCascade.cxx:582
IVertexFitter.h
EventInfo.h
MuonContainer.h
VxCascadeInfo.h
mc.mass_b
mass_b
Definition: mc.PhPy8EG_A14NNPDF23_gg4l_example.py:21
DerivationFramework::VertexLinkVector
std::vector< VertexLink > VertexLinkVector
Definition: Cascade3Plus1.cxx:24
VertexContainer.h
xAOD::Vertex_v1
Class describing a Vertex.
Definition: Vertex_v1.h:42
DerivationFramework::MuPlusDsCascade::m_iVertexFitter
ToolHandle< Trk::TrkVKalVrtFitter > m_iVertexFitter
Definition: MuPlusDsCascade.h:78
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
DerivationFramework::MuPlusDsCascade::m_MassUpper
double m_MassUpper
Definition: MuPlusDsCascade.h:65
xAOD::Vertex_v1::y
float y() const
Returns the y position.
DerivationFramework::MuPlusDsCascade::addBranches
virtual StatusCode addBranches() const override
Pass the thinning service
Definition: MuPlusDsCascade.cxx:95
SG::ConstAccessor< T, AuxAllocator_t< T > >::isAvailable
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
DerivationFramework::MuPlusDsCascade::m_refPVContainerName
std::string m_refPVContainerName
Definition: MuPlusDsCascade.h:86
xAOD::TrackParticle_v1
Class describing a TrackParticle.
Definition: TrackParticle_v1.h:43
AthAlgTool
Definition: AthAlgTool.h:26
ITrackSelectorTool.h
CaloNoise_fillDB.mu
mu
Definition: CaloNoise_fillDB.py:53
DerivationFramework::MuPlusDsCascade::m_vtx1Daug1MassHypo
double m_vtx1Daug1MassHypo
Definition: MuPlusDsCascade.h:69
DataVector::size
size_type size() const noexcept
Returns the number of elements in the collection.
DerivationFramework::MuPlusDsCascade::m_hypoName
std::string m_hypoName
name of the mass hypothesis.
Definition: MuPlusDsCascade.h:87
DerivationFramework::MuPlusDsCascade::m_Dx_pid
int m_Dx_pid
Definition: MuPlusDsCascade.h:73
DerivationFramework::MuPlusDsCascade::m_DxMassLower
double m_DxMassLower
Definition: MuPlusDsCascade.h:62
HepMCHelpers.h
xAOD::BPhysHypoHelper::setMassErr
bool setMassErr(const float val)
invariant mass error
Definition: BPhysHypoHelper.cxx:54
VertexAuxContainer.h
xAOD::TrackParticle_v1::phi
virtual double phi() const override final
The azimuthal angle ( ) of the particle (has range to .)