ATLAS Offline Software
JpsiPlusDpstCascade.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
5 // JpsiPlusDpstCascade.cxx, (c) ATLAS Detector software
7 #include "DerivationFrameworkBPhys/JpsiPlusDpstCascade.h"
8 #include "TrkVertexFitterInterfaces/IVertexFitter.h"
9 #include "TrkVKalVrtFitter/TrkVKalVrtFitter.h"
10 #include "TrkVertexAnalysisUtils/V0Tools.h"
11 #include "GaudiKernel/IPartPropSvc.h"
12 #include "DerivationFrameworkBPhys/CascadeTools.h"
13 #include "DerivationFrameworkBPhys/BPhysPVCascadeTools.h"
14 #include "xAODTracking/VertexAuxContainer.h"
15 #include "xAODBPhys/BPhysHypoHelper.h"
16 #include <algorithm>
17 #include "xAODTracking/VertexContainer.h"
18 #include "HepPDT/ParticleDataTable.hh"
19 #include "TruthUtils/HepMCHelpers.h"
20 
21 #include "TrkVKalVrtFitter/VxCascadeInfo.h"
22 
23 namespace DerivationFramework {
24  typedef ElementLink<xAOD::VertexContainer> VertexLink;
25  typedef std::vector<VertexLink> VertexLinkVector;
26  typedef std::vector<const xAOD::TrackParticle*> TrackBag;
27 
28 
29  StatusCode JpsiPlusDpstCascade::initialize() {
30 
31  // retrieving vertex Fitter
32  ATH_CHECK( m_iVertexFitter.retrieve());
33 
34  // retrieving the V0 tools
35  ATH_CHECK( m_V0Tools.retrieve());
36 
37  // retrieving the Cascade tools
38  ATH_CHECK( m_CascadeTools.retrieve());
39 
40  // Get the beam spot service
41  ATH_CHECK(m_eventInfo_key.initialize());
42 
43  IPartPropSvc* partPropSvc = nullptr;
44  ATH_CHECK( service("PartPropSvc", partPropSvc, true) );
45  auto pdt = partPropSvc->PDT();
46 
47  // retrieve particle masses
48  if(m_mass_jpsi < 0. ) m_mass_jpsi = BPhysPVCascadeTools::getParticleMass(pdt, MC::JPSI);
49  if(m_vtx0MassHypo < 0.) m_vtx0MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::BCPLUS);
50  if(m_vtx1MassHypo < 0.) m_vtx1MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::D0);
51 
52  if(m_vtx0Daug1MassHypo < 0.) m_vtx0Daug1MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::MUON);
53  if(m_vtx0Daug2MassHypo < 0.) m_vtx0Daug2MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::MUON);
54  if(m_vtx0Daug3MassHypo < 0.) m_vtx0Daug3MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::PIPLUS);
55  if(m_vtx1Daug1MassHypo < 0.) m_vtx1Daug1MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::PIPLUS);
56  if(m_vtx1Daug2MassHypo < 0.) m_vtx1Daug2MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::KPLUS);
57 
58  return StatusCode::SUCCESS;
59  }
60 
61 
62  StatusCode JpsiPlusDpstCascade::addBranches() const
63  {
64  std::vector<Trk::VxCascadeInfo*> cascadeinfoContainer;
65  constexpr int topoN = 2;
66  std::array<xAOD::VertexContainer*, topoN> Vtxwritehandles;
67  std::array<xAOD::VertexAuxContainer*, topoN> Vtxwritehandlesaux;
68  if(m_cascadeOutputsKeys.size() !=topoN) { ATH_MSG_FATAL("Incorrect number of VtxContainers"); return StatusCode::FAILURE; }
69 
70  for(int i =0; i<topoN;i++){
71  Vtxwritehandles[i] = new xAOD::VertexContainer();
72  Vtxwritehandlesaux[i] = new xAOD::VertexAuxContainer();
73  Vtxwritehandles[i]->setStore(Vtxwritehandlesaux[i]);
74  ATH_CHECK(evtStore()->record(Vtxwritehandles[i] , m_cascadeOutputsKeys[i] ));
75  ATH_CHECK(evtStore()->record(Vtxwritehandlesaux[i], m_cascadeOutputsKeys[i] + "Aux."));
76  }
77 
78  //----------------------------------------------------
79  // retrieve primary vertices
80  //----------------------------------------------------
81  const xAOD::Vertex * primaryVertex(nullptr);
82  const xAOD::VertexContainer *pvContainer(nullptr);
83  ATH_CHECK(evtStore()->retrieve(pvContainer, m_VxPrimaryCandidateName));
84  ATH_MSG_DEBUG("Found " << m_VxPrimaryCandidateName << " in StoreGate!");
85 
86  if (pvContainer->size()==0){
87  ATH_MSG_WARNING("You have no primary vertices: " << pvContainer->size());
88  return StatusCode::RECOVERABLE;
89  } else {
90  primaryVertex = (*pvContainer)[0];
91  }
92 
93  //----------------------------------------------------
94  // Try to retrieve refitted primary vertices
95  //----------------------------------------------------
96  xAOD::VertexContainer* refPvContainer = nullptr;
97  xAOD::VertexAuxContainer* refPvAuxContainer = nullptr;
98  if (m_refitPV) {
99  if (evtStore()->contains<xAOD::VertexContainer>(m_refPVContainerName)) {
100  // refitted PV container exists. Get it from the store gate
101  ATH_CHECK(evtStore()->retrieve(refPvContainer , m_refPVContainerName ));
102  ATH_CHECK(evtStore()->retrieve(refPvAuxContainer, m_refPVContainerName + "Aux."));
103  } else {
104  // refitted PV container does not exist. Create a new one.
105  refPvContainer = new xAOD::VertexContainer;
106  refPvAuxContainer = new xAOD::VertexAuxContainer;
107  refPvContainer->setStore(refPvAuxContainer);
108  ATH_CHECK(evtStore()->record(refPvContainer , m_refPVContainerName));
109  ATH_CHECK(evtStore()->record(refPvAuxContainer, m_refPVContainerName+"Aux."));
110  }
111  }
112 
113  ATH_CHECK(performSearch(&cascadeinfoContainer));
114 
115  SG::ReadHandle<xAOD::EventInfo> evt(m_eventInfo_key);
116  if(not evt.isValid()) ATH_MSG_ERROR("Cannot Retrieve " << evt.key() );
117  BPhysPVCascadeTools helper(&(*m_CascadeTools), evt.cptr());
118  helper.SetMinNTracksInPV(m_PV_minNTracks);
119 
120  // Decorators for the main vertex: chi2, ndf, pt and pt error, plus the D0 vertex variables
121  SG::AuxElement::Decorator<VertexLinkVector> CascadeLinksDecor("CascadeVertexLinks");
122  SG::AuxElement::Decorator<VertexLinkVector> JpsipiLinksDecor("JpsipiVertexLinks");
123  SG::AuxElement::Decorator<VertexLinkVector> D0LinksDecor("D0VertexLinks");
124  SG::AuxElement::Decorator<float> chi2_decor("ChiSquared");
125  SG::AuxElement::Decorator<float> ndof_decor("NumberDoF");
126  SG::AuxElement::Decorator<float> Pt_decor("Pt");
127  SG::AuxElement::Decorator<float> PtErr_decor("PtErr");
128  SG::AuxElement::Decorator<float> Mass_svdecor("D0_mass");
129  SG::AuxElement::Decorator<float> MassErr_svdecor("D0_massErr");
130  SG::AuxElement::Decorator<float> Pt_svdecor("D0_Pt");
131  SG::AuxElement::Decorator<float> PtErr_svdecor("D0_PtErr");
132  SG::AuxElement::Decorator<float> Lxy_svdecor("D0_Lxy");
133  SG::AuxElement::Decorator<float> LxyErr_svdecor("D0_LxyErr");
134  SG::AuxElement::Decorator<float> Tau_svdecor("D0_Tau");
135  SG::AuxElement::Decorator<float> TauErr_svdecor("D0_TauErr");
136 
137  SG::AuxElement::Decorator<float> MassMumu_decor("Mumu_mass");
138  SG::AuxElement::Decorator<float> MassKpi_svdecor("Kpi_mass");
139  SG::AuxElement::Decorator<float> MassJpsi_decor("Jpsi_mass");
140  SG::AuxElement::Decorator<float> MassPiD0_decor("PiD0_mass");
141 
142  ATH_MSG_DEBUG("cascadeinfoContainer size " << cascadeinfoContainer.size());
143 
144  // Get Jpsi+pi container and identify the input Jpsi+pi
145  const xAOD::VertexContainer *jpsipiContainer(nullptr);
146  ATH_CHECK(evtStore()->retrieve(jpsipiContainer , m_vertexContainerKey ));
147  const xAOD::VertexContainer *d0Container(nullptr);
148  ATH_CHECK(evtStore()->retrieve(d0Container , m_vertexD0ContainerKey ));
149 
150  for (Trk::VxCascadeInfo* x : cascadeinfoContainer) {
151  if(x==nullptr) {
152  ATH_MSG_ERROR("cascadeinfoContainer is null");
153  //x is dereferenced if we pass this
154  return StatusCode::FAILURE;
155  }
156 
157  // the cascade fitter returns:
158  // std::vector<xAOD::Vertex*>, each xAOD::Vertex contains the refitted track parameters (perigee at the vertex position)
159  // vertices[iv] the links to the original TPs and a covariance of size 3+5*NTRK; the chi2 of the total fit
160  // is split between the cascade vertices as per track contribution
161  // std::vector< std::vector<TLorentzVector> >, each std::vector<TLorentzVector> contains the refitted momenta (TLorentzVector)
162  // momenta[iv][...] of all tracks in the corresponding vertex, including any pseudotracks (from cascade vertices)
163  // originating in this vertex; the masses are as assigned in the cascade fit
164  // std::vector<Amg::MatrixX>, the corresponding covariance matrices in momentum space
165  // covariance[iv]
166  // int nDoF, double Chi2
167  //
168  // the invariant mass, pt, lifetime etc. errors should be calculated using the covariance matrices in momentum space as these
169  // take into account the full track-track and track-vertex correlations
170  //
171  // 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.
172  // The covariance terms between the two vertices are not stored. In momentum space momenta[0] contains the 2 V0 tracks,
173  // their momenta add up to the momentum of the 3rd track in momenta[1], the first two being the Jpsi tracks
174 
175  const std::vector<xAOD::Vertex*> &cascadeVertices = x->vertices();
176  if(cascadeVertices.size()!=topoN)
177  ATH_MSG_ERROR("Incorrect number of vertices");
178  if(cascadeVertices[0] == nullptr || cascadeVertices[1] == nullptr) ATH_MSG_ERROR("Error null vertex");
179  // Keep vertices (bear in mind that they come in reverse order!)
180  for(int i =0;i<topoN;i++) Vtxwritehandles[i]->push_back(cascadeVertices[i]);
181 
182  x->setSVOwnership(false); // Prevent Container from deleting vertices
183  const auto mainVertex = cascadeVertices[1]; // this is the B_c+/- vertex
184  const std::vector< std::vector<TLorentzVector> > &moms = x->getParticleMoms();
185 
186  // Set links to cascade vertices
187  std::vector<const xAOD::Vertex*> verticestoLink;
188  verticestoLink.push_back(cascadeVertices[0]);
189  if(Vtxwritehandles[1] == nullptr) ATH_MSG_ERROR("Vtxwritehandles[1] is null");
190  if(!BPhysPVCascadeTools::LinkVertices(CascadeLinksDecor, verticestoLink, Vtxwritehandles[0], cascadeVertices[1]))
191  ATH_MSG_ERROR("Error decorating with cascade vertices");
192 
193  // Identify the input Jpsi+pi
194  const xAOD::Vertex* jpsipiVertex = BPhysPVCascadeTools::FindVertex<3>(jpsipiContainer, cascadeVertices[1]);
195  ATH_MSG_DEBUG("1 pt Jpsi+pi tracks " << cascadeVertices[1]->trackParticle(0)->pt() << ", " << cascadeVertices[1]->trackParticle(1)->pt() << ", " << cascadeVertices[1]->trackParticle(2)->pt());
196  if (jpsipiVertex) ATH_MSG_DEBUG("2 pt Jpsi+pi tracks " << jpsipiVertex->trackParticle(0)->pt() << ", " << jpsipiVertex->trackParticle(1)->pt() << ", " << jpsipiVertex->trackParticle(2)->pt());
197 
198  // Identify the input D0
199  const xAOD::Vertex* d0Vertex = BPhysPVCascadeTools::FindVertex<2>(d0Container, cascadeVertices[0]);;
200  ATH_MSG_DEBUG("1 pt D0 tracks " << cascadeVertices[0]->trackParticle(0)->pt() << ", " << cascadeVertices[0]->trackParticle(1)->pt());
201  if (d0Vertex) ATH_MSG_DEBUG("2 pt D0 tracks " << d0Vertex->trackParticle(0)->pt() << ", " << d0Vertex->trackParticle(1)->pt());
202 
203  // Set links to input vertices
204  std::vector<const xAOD::Vertex*> jpsipiVerticestoLink;
205  if (jpsipiVertex) jpsipiVerticestoLink.push_back(jpsipiVertex);
206  else ATH_MSG_WARNING("Could not find linking Jpsi+pi");
207  if(!BPhysPVCascadeTools::LinkVertices(JpsipiLinksDecor, jpsipiVerticestoLink, jpsipiContainer, cascadeVertices[1]))
208  ATH_MSG_ERROR("Error decorating with Jpsi+pi vertices");
209 
210  std::vector<const xAOD::Vertex*> d0VerticestoLink;
211  if (d0Vertex) d0VerticestoLink.push_back(d0Vertex);
212  else ATH_MSG_WARNING("Could not find linking D0");
213  if(!BPhysPVCascadeTools::LinkVertices(D0LinksDecor, d0VerticestoLink, d0Container, cascadeVertices[1]))
214  ATH_MSG_ERROR("Error decorating with D0 vertices");
215 
216  bool tagD0(true);
217  if (jpsipiVertex){
218  if(abs(m_Dx_pid)==421 && (jpsipiVertex->trackParticle(2)->charge()==-1)) tagD0 = false;
219  }
220 
221  double mass_b = m_vtx0MassHypo;
222  double mass_d0 = m_vtx1MassHypo;
223  std::vector<double> massesJpsipi;
224  massesJpsipi.push_back(m_vtx0Daug1MassHypo);
225  massesJpsipi.push_back(m_vtx0Daug2MassHypo);
226  massesJpsipi.push_back(m_vtx0Daug3MassHypo);
227  std::vector<double> massesD0;
228  if(tagD0){
229  massesD0.push_back(m_vtx1Daug1MassHypo);
230  massesD0.push_back(m_vtx1Daug2MassHypo);
231  }else{ // Change the oreder of masses for D*-->D0bar pi-, D0bar->K+pi-
232  massesD0.push_back(m_vtx1Daug2MassHypo);
233  massesD0.push_back(m_vtx1Daug1MassHypo);
234  }
235  std::vector<double> Masses;
236  Masses.push_back(m_vtx0Daug1MassHypo);
237  Masses.push_back(m_vtx0Daug2MassHypo);
238  Masses.push_back(m_vtx0Daug3MassHypo);
239  Masses.push_back(m_vtx1MassHypo);
240 
241  // loop over candidates -- Don't apply PV_minNTracks requirement here
242  // because it may result in exclusion of the high-pt PV.
243  // get good PVs
244 
245  xAOD::BPhysHypoHelper vtx(m_hypoName, mainVertex);
246 
247  // Get refitted track momenta from all vertices, charged tracks only
248  BPhysPVCascadeTools::SetVectorInfo(vtx, x);
249 
250  // Decorate main vertex
251  //
252  // 1.a) mass, mass error
253  BPHYS_CHECK( vtx.setMass(m_CascadeTools->invariantMass(moms[1])) );
254  BPHYS_CHECK( vtx.setMassErr(m_CascadeTools->invariantMassError(moms[1],x->getCovariance()[1])) );
255  // 1.b) pt and pT error (the default pt of mainVertex is != the pt of the full cascade fit!)
256  Pt_decor(*mainVertex) = m_CascadeTools->pT(moms[1]);
257  PtErr_decor(*mainVertex) = m_CascadeTools->pTError(moms[1],x->getCovariance()[1]);
258  // 1.c) chi2 and ndof (the default chi2 of mainVertex is != the chi2 of the full cascade fit!)
259  chi2_decor(*mainVertex) = x->fitChi2();
260  ndof_decor(*mainVertex) = x->nDoF();
261 
262  float massMumu = 0.;
263  if (jpsipiVertex) {
264  TLorentzVector p4_mu1, p4_mu2;
265  p4_mu1.SetPtEtaPhiM(jpsipiVertex->trackParticle(0)->pt(),
266  jpsipiVertex->trackParticle(0)->eta(),
267  jpsipiVertex->trackParticle(0)->phi(), m_vtx0Daug1MassHypo);
268  p4_mu2.SetPtEtaPhiM(jpsipiVertex->trackParticle(1)->pt(),
269  jpsipiVertex->trackParticle(1)->eta(),
270  jpsipiVertex->trackParticle(1)->phi(), m_vtx0Daug2MassHypo);
271  massMumu = (p4_mu1 + p4_mu2).M();
272  }
273  MassMumu_decor(*mainVertex) = massMumu;
274 
275  float massKpi = 0.;
276  if (d0Vertex) {
277  TLorentzVector p4_ka, p4_pi;
278  if(tagD0){
279  p4_pi.SetPtEtaPhiM(d0Vertex->trackParticle(0)->pt(),
280  d0Vertex->trackParticle(0)->eta(),
281  d0Vertex->trackParticle(0)->phi(), m_vtx1Daug1MassHypo);
282  p4_ka.SetPtEtaPhiM(d0Vertex->trackParticle(1)->pt(),
283  d0Vertex->trackParticle(1)->eta(),
284  d0Vertex->trackParticle(1)->phi(), m_vtx1Daug2MassHypo);
285  }else{ // Change the oreder of masses for D*-->D0bar pi-, D0bar->K+pi-
286  p4_pi.SetPtEtaPhiM(d0Vertex->trackParticle(1)->pt(),
287  d0Vertex->trackParticle(1)->eta(),
288  d0Vertex->trackParticle(1)->phi(), m_vtx1Daug1MassHypo);
289  p4_ka.SetPtEtaPhiM(d0Vertex->trackParticle(0)->pt(),
290  d0Vertex->trackParticle(0)->eta(),
291  d0Vertex->trackParticle(0)->phi(), m_vtx1Daug2MassHypo);
292  }
293  massKpi = (p4_ka + p4_pi).M();
294  }
295  MassKpi_svdecor(*mainVertex) = massKpi;
296  MassJpsi_decor(*mainVertex) = (moms[1][0] + moms[1][1]).M();
297  MassPiD0_decor(*mainVertex) = (moms[1][2] + moms[1][3]).M();
298 
299 
300  ATH_CHECK(helper.FillCandwithRefittedVertices(m_refitPV, pvContainer,
301  refPvContainer, &(*m_pvRefitter), m_PV_max, m_DoVertexType, x, 1, mass_b, vtx));
302 
303  // 4) decorate the main vertex with D0 vertex mass, pt, lifetime and lxy values (plus errors)
304  // D0 points to the main vertex, so lifetime and lxy are w.r.t the main vertex
305  Mass_svdecor(*mainVertex) = m_CascadeTools->invariantMass(moms[0]);
306  MassErr_svdecor(*mainVertex) = m_CascadeTools->invariantMassError(moms[0],x->getCovariance()[0]);
307  Pt_svdecor(*mainVertex) = m_CascadeTools->pT(moms[0]);
308  PtErr_svdecor(*mainVertex) = m_CascadeTools->pTError(moms[0],x->getCovariance()[0]);
309  Lxy_svdecor(*mainVertex) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
310  LxyErr_svdecor(*mainVertex) = m_CascadeTools->lxyError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
311  Tau_svdecor(*mainVertex) = m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1]);
312  TauErr_svdecor(*mainVertex) = m_CascadeTools->tauError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
313 
314  // Some checks in DEBUG mode
315  ATH_MSG_DEBUG("chi2 " << x->fitChi2()
316  << " chi2_1 " << m_V0Tools->chisq(cascadeVertices[0])
317  << " chi2_2 " << m_V0Tools->chisq(cascadeVertices[1])
318  << " vprob " << m_CascadeTools->vertexProbability(x->nDoF(),x->fitChi2()));
319  ATH_MSG_DEBUG("ndf " << x->nDoF() << " ndf_1 " << m_V0Tools->ndof(cascadeVertices[0]) << " ndf_2 " << m_V0Tools->ndof(cascadeVertices[1]));
320  ATH_MSG_DEBUG("V0Tools mass_d0 " << m_V0Tools->invariantMass(cascadeVertices[0],massesD0)
321  << " error " << m_V0Tools->invariantMassError(cascadeVertices[0],massesD0)
322  << " mass_J " << m_V0Tools->invariantMass(cascadeVertices[1],massesJpsipi)
323  << " error " << m_V0Tools->invariantMassError(cascadeVertices[1],massesJpsipi));
324  // masses and errors, using track masses assigned in the fit
325  double Mass_B = m_CascadeTools->invariantMass(moms[1]);
326  double Mass_D0 = m_CascadeTools->invariantMass(moms[0]);
327  double Mass_B_err = m_CascadeTools->invariantMassError(moms[1],x->getCovariance()[1]);
328  double Mass_D0_err = m_CascadeTools->invariantMassError(moms[0],x->getCovariance()[0]);
329  ATH_MSG_DEBUG("Mass_B " << Mass_B << " Mass_D0 " << Mass_D0);
330  ATH_MSG_DEBUG("Mass_B_err " << Mass_B_err << " Mass_D0_err " << Mass_D0_err);
331  double mprob_B = m_CascadeTools->massProbability(mass_b,Mass_B,Mass_B_err);
332  double mprob_D0 = m_CascadeTools->massProbability(mass_d0,Mass_D0,Mass_D0_err);
333  ATH_MSG_DEBUG("mprob_B " << mprob_B << " mprob_D0 " << mprob_D0);
334  // masses and errors, assigning user defined track masses
335  ATH_MSG_DEBUG("Mass_b " << m_CascadeTools->invariantMass(moms[1],Masses)
336  << " Mass_d0 " << m_CascadeTools->invariantMass(moms[0],massesD0));
337  ATH_MSG_DEBUG("Mass_b_err " << m_CascadeTools->invariantMassError(moms[1],x->getCovariance()[1],Masses)
338  << " Mass_d0_err " << m_CascadeTools->invariantMassError(moms[0],x->getCovariance()[0],massesD0));
339  ATH_MSG_DEBUG("pt_b " << m_CascadeTools->pT(moms[1])
340  << " pt_d " << m_CascadeTools->pT(moms[0])
341  << " pt_d0 " << m_V0Tools->pT(cascadeVertices[0]));
342  ATH_MSG_DEBUG("ptErr_b " << m_CascadeTools->pTError(moms[1],x->getCovariance()[1])
343  << " ptErr_d " << m_CascadeTools->pTError(moms[0],x->getCovariance()[0])
344  << " ptErr_d0 " << m_V0Tools->pTError(cascadeVertices[0]));
345  ATH_MSG_DEBUG("lxy_B " << m_V0Tools->lxy(cascadeVertices[1],primaryVertex) << " lxy_D " << m_V0Tools->lxy(cascadeVertices[0],cascadeVertices[1]));
346  ATH_MSG_DEBUG("lxy_b " << m_CascadeTools->lxy(moms[1],cascadeVertices[1],primaryVertex) << " lxy_d " << m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]));
347  ATH_MSG_DEBUG("lxyErr_b " << m_CascadeTools->lxyError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
348  << " lxyErr_d " << m_CascadeTools->lxyError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
349  << " lxyErr_d0 " << m_V0Tools->lxyError(cascadeVertices[0],cascadeVertices[1]));
350  ATH_MSG_DEBUG("tau_B " << m_CascadeTools->tau(moms[1],cascadeVertices[1],primaryVertex,mass_b)
351  << " tau_d0 " << m_V0Tools->tau(cascadeVertices[0],cascadeVertices[1],massesD0));
352  ATH_MSG_DEBUG("tau_b " << m_CascadeTools->tau(moms[1],cascadeVertices[1],primaryVertex)
353  << " tau_d " << m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1])
354  << " tau_D " << m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1],mass_d0));
355  ATH_MSG_DEBUG("tauErr_b " << m_CascadeTools->tauError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
356  << " tauErr_d " << m_CascadeTools->tauError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
357  << " tauErr_d0 " << m_V0Tools->tauError(cascadeVertices[0],cascadeVertices[1],massesD0));
358  ATH_MSG_DEBUG("TauErr_b " << m_CascadeTools->tauError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex,mass_b)
359  << " TauErr_d " << m_CascadeTools->tauError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1],mass_d0)
360  << " TauErr_d0 " << m_V0Tools->tauError(cascadeVertices[0],cascadeVertices[1],massesD0,mass_d0));
361 
362  ATH_MSG_DEBUG("CascadeTools main vert wrt PV " << " CascadeTools SV " << " V0Tools SV");
363  ATH_MSG_DEBUG("a0z " << m_CascadeTools->a0z(moms[1],cascadeVertices[1],primaryVertex)
364  << ", " << m_CascadeTools->a0z(moms[0],cascadeVertices[0],cascadeVertices[1])
365  << ", " << m_V0Tools->a0z(cascadeVertices[0],cascadeVertices[1]));
366  ATH_MSG_DEBUG("a0zErr " << m_CascadeTools->a0zError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
367  << ", " << m_CascadeTools->a0zError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
368  << ", " << m_V0Tools->a0zError(cascadeVertices[0],cascadeVertices[1]));
369  ATH_MSG_DEBUG("a0xy " << m_CascadeTools->a0xy(moms[1],cascadeVertices[1],primaryVertex)
370  << ", " << m_CascadeTools->a0xy(moms[0],cascadeVertices[0],cascadeVertices[1])
371  << ", " << m_V0Tools->a0xy(cascadeVertices[0],cascadeVertices[1]));
372  ATH_MSG_DEBUG("a0xyErr " << m_CascadeTools->a0xyError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
373  << ", " << m_CascadeTools->a0xyError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
374  << ", " << m_V0Tools->a0xyError(cascadeVertices[0],cascadeVertices[1]));
375  ATH_MSG_DEBUG("a0 " << m_CascadeTools->a0(moms[1],cascadeVertices[1],primaryVertex)
376  << ", " << m_CascadeTools->a0(moms[0],cascadeVertices[0],cascadeVertices[1])
377  << ", " << m_V0Tools->a0(cascadeVertices[0],cascadeVertices[1]));
378  ATH_MSG_DEBUG("a0Err " << m_CascadeTools->a0Error(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
379  << ", " << m_CascadeTools->a0Error(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
380  << ", " << m_V0Tools->a0Error(cascadeVertices[0],cascadeVertices[1]));
381  ATH_MSG_DEBUG("x0 " << m_V0Tools->vtx(cascadeVertices[0]).x() << " y0 " << m_V0Tools->vtx(cascadeVertices[0]).y() << " z0 " << m_V0Tools->vtx(cascadeVertices[0]).z());
382  ATH_MSG_DEBUG("x1 " << m_V0Tools->vtx(cascadeVertices[1]).x() << " y1 " << m_V0Tools->vtx(cascadeVertices[1]).y() << " z1 " << m_V0Tools->vtx(cascadeVertices[1]).z());
383  ATH_MSG_DEBUG("X0 " << primaryVertex->x() << " Y0 " << primaryVertex->y() << " Z0 " << primaryVertex->z());
384  ATH_MSG_DEBUG("rxy0 " << m_V0Tools->rxy(cascadeVertices[0]) << " rxyErr0 " << m_V0Tools->rxyError(cascadeVertices[0]));
385  ATH_MSG_DEBUG("rxy1 " << m_V0Tools->rxy(cascadeVertices[1]) << " rxyErr1 " << m_V0Tools->rxyError(cascadeVertices[1]));
386  ATH_MSG_DEBUG("Rxy0 wrt PV " << m_V0Tools->rxy(cascadeVertices[0],primaryVertex) << " RxyErr0 wrt PV " << m_V0Tools->rxyError(cascadeVertices[0],primaryVertex));
387  ATH_MSG_DEBUG("Rxy1 wrt PV " << m_V0Tools->rxy(cascadeVertices[1],primaryVertex) << " RxyErr1 wrt PV " << m_V0Tools->rxyError(cascadeVertices[1],primaryVertex));
388  ATH_MSG_DEBUG("number of covariance matrices " << (x->getCovariance()).size());
389  } // loop over cascadeinfoContainer
390 
391  // Deleting cascadeinfo since this won't be stored.
392  // Vertices have been kept in m_cascadeOutputs and should be owned by their container
393  for (auto x : cascadeinfoContainer) delete x;
394 
395  return StatusCode::SUCCESS;
396  }
397 
398 
399  JpsiPlusDpstCascade::JpsiPlusDpstCascade(const std::string& t, const std::string& n, const IInterface* p) : AthAlgTool(t,n,p),
400  m_vertexContainerKey(""),
401  m_vertexD0ContainerKey(""),
402  m_cascadeOutputsKeys{ "JpsiPlusDpstCascadeVtx1", "JpsiPlusDpstCascadeVtx2" },
403  m_VxPrimaryCandidateName("PrimaryVertices"),
404  m_jpsiMassLower(0.0),
405  m_jpsiMassUpper(10000.0),
406  m_jpsipiMassLower(0.0),
407  m_jpsipiMassUpper(10000.0),
408  m_D0MassLower(0.0),
409  m_D0MassUpper(10000.0),
410  m_DstMassLower(0.0),
411  m_DstMassUpper(10000.0),
412  m_MassLower(0.0),
413  m_MassUpper(20000.0),
414  m_vtx0MassHypo(-1),
415  m_vtx1MassHypo(-1),
416  m_vtx0Daug1MassHypo(-1),
417  m_vtx0Daug2MassHypo(-1),
418  m_vtx0Daug3MassHypo(-1),
419  m_vtx1Daug1MassHypo(-1),
420  m_vtx1Daug2MassHypo(-1),
421  m_mass_jpsi(-1),
422  m_Dx_pid(421),
423  m_constrD0(true),
424  m_constrJpsi(true),
425  m_chi2cut(-1.0),
426  m_iVertexFitter("Trk::TrkVKalVrtFitter"),
427  m_pvRefitter("Analysis::PrimaryVertexRefitter"),
428  m_V0Tools("Trk::V0Tools"),
429  m_CascadeTools("DerivationFramework::CascadeTools")
430  {
431  declareProperty("JpsipiVertices", m_vertexContainerKey);
432  declareProperty("D0Vertices", m_vertexD0ContainerKey);
433  declareProperty("VxPrimaryCandidateName", m_VxPrimaryCandidateName);
434  declareProperty("RefPVContainerName", m_refPVContainerName = "RefittedPrimaryVertices");
435  declareProperty("JpsiMassLowerCut", m_jpsiMassLower);
436  declareProperty("JpsiMassUpperCut", m_jpsiMassUpper);
437  declareProperty("JpsipiMassLowerCut", m_jpsipiMassLower);
438  declareProperty("JpsipiMassUpperCut", m_jpsipiMassUpper);
439  declareProperty("D0MassLowerCut", m_D0MassLower);
440  declareProperty("D0MassUpperCut", m_D0MassUpper);
441  declareProperty("DstMassLowerCut", m_DstMassLower);
442  declareProperty("DstMassUpperCut", m_DstMassUpper);
443  declareProperty("MassLowerCut", m_MassLower);
444  declareProperty("MassUpperCut", m_MassUpper);
445  declareProperty("HypothesisName", m_hypoName = "Bc");
446  declareProperty("Vtx0MassHypo", m_vtx0MassHypo);
447  declareProperty("Vtx1MassHypo", m_vtx1MassHypo);
448  declareProperty("Vtx0Daug1MassHypo", m_vtx0Daug1MassHypo);
449  declareProperty("Vtx0Daug2MassHypo", m_vtx0Daug2MassHypo);
450  declareProperty("Vtx0Daug3MassHypo", m_vtx0Daug3MassHypo);
451  declareProperty("Vtx1Daug1MassHypo", m_vtx1Daug1MassHypo);
452  declareProperty("Vtx1Daug2MassHypo", m_vtx1Daug2MassHypo);
453  declareProperty("JpsiMass", m_mass_jpsi);
454  declareProperty("DxHypothesis", m_Dx_pid);
455  declareProperty("ApplyD0MassConstraint", m_constrD0);
456  declareProperty("ApplyJpsiMassConstraint", m_constrJpsi);
457  declareProperty("Chi2Cut", m_chi2cut);
458  declareProperty("RefitPV", m_refitPV = true);
459  declareProperty("MaxnPV", m_PV_max = 999);
460  declareProperty("MinNTracksInPV", m_PV_minNTracks = 0);
461  declareProperty("DoVertexType", m_DoVertexType = 7);
462  declareProperty("TrkVertexFitterTool", m_iVertexFitter);
463  declareProperty("PVRefitter", m_pvRefitter);
464  declareProperty("V0Tools", m_V0Tools);
465  declareProperty("CascadeTools", m_CascadeTools);
466  declareProperty("CascadeVertexCollections", m_cascadeOutputsKeys);
467  }
468 
469  JpsiPlusDpstCascade::~JpsiPlusDpstCascade(){ }
470 
471  StatusCode JpsiPlusDpstCascade::performSearch(std::vector<Trk::VxCascadeInfo*> *cascadeinfoContainer) const
472  {
473  ATH_MSG_DEBUG( "JpsiPlusDpstCascade::performSearch" );
474  assert(cascadeinfoContainer!=nullptr);
475 
476  // Get TrackParticle container (for setting links to the original tracks)
477  const xAOD::TrackParticleContainer *trackContainer(nullptr);
478  ATH_CHECK(evtStore()->retrieve(trackContainer , "InDetTrackParticles" ));
479 
480  // Get Jpsi+pi container
481  const xAOD::VertexContainer *jpsipiContainer(nullptr);
482  ATH_CHECK(evtStore()->retrieve(jpsipiContainer , m_vertexContainerKey ));
483 
484  // Get D0 container
485  const xAOD::VertexContainer *d0Container(nullptr);
486  ATH_CHECK(evtStore()->retrieve(d0Container , m_vertexD0ContainerKey ));
487 
488  double mass_d0 = m_vtx1MassHypo;
489  std::vector<const xAOD::TrackParticle*> tracksJpsipi;
490  std::vector<const xAOD::TrackParticle*> tracksJpsi;
491  std::vector<const xAOD::TrackParticle*> tracksD0;
492  std::vector<const xAOD::TrackParticle*> tracksBc;
493  std::vector<double> massesJpsipi;
494  massesJpsipi.push_back(m_vtx0Daug1MassHypo);
495  massesJpsipi.push_back(m_vtx0Daug2MassHypo);
496  massesJpsipi.push_back(m_vtx0Daug3MassHypo);
497  std::vector<double> massesD0;
498  massesD0.push_back(m_vtx1Daug1MassHypo);
499  massesD0.push_back(m_vtx1Daug2MassHypo);
500  std::vector<double> massesD0b; // Change the oreder of masses for D*-->D0bar pi-, D0bar->K+pi-
501  massesD0b.push_back(m_vtx1Daug2MassHypo);
502  massesD0b.push_back(m_vtx1Daug1MassHypo);
503  std::vector<double> Masses;
504  Masses.push_back(m_vtx0Daug1MassHypo);
505  Masses.push_back(m_vtx0Daug2MassHypo);
506  Masses.push_back(m_vtx0Daug3MassHypo);
507  Masses.push_back(m_vtx1MassHypo);
508 
509  // Select J/psi pi+ candidates before calling cascade fit
510  std::vector<const xAOD::Vertex*> selectedJpsipiCandidates;
511  for(auto vxcItr=jpsipiContainer->cbegin(); vxcItr!=jpsipiContainer->cend(); ++vxcItr) {
512 
513  // Check the passed flag first
514  const xAOD::Vertex* vtx = *vxcItr;
515  SG::AuxElement::Accessor<Char_t> flagAcc1("passed_Jpsipi");
516  if(flagAcc1.isAvailable(*vtx)){
517  if(!flagAcc1(*vtx)) continue;
518  }
519 
520  // Check J/psi candidate invariant mass and skip if need be
521  TLorentzVector p4Mup_in, p4Mum_in;
522  p4Mup_in.SetPtEtaPhiM((*vxcItr)->trackParticle(0)->pt(),
523  (*vxcItr)->trackParticle(0)->eta(),
524  (*vxcItr)->trackParticle(0)->phi(), m_vtx0Daug1MassHypo);
525  p4Mum_in.SetPtEtaPhiM((*vxcItr)->trackParticle(1)->pt(),
526  (*vxcItr)->trackParticle(1)->eta(),
527  (*vxcItr)->trackParticle(1)->phi(), m_vtx0Daug2MassHypo);
528  double mass_Jpsi = (p4Mup_in + p4Mum_in).M();
529  ATH_MSG_DEBUG("Jpsi mass " << mass_Jpsi);
530  if (mass_Jpsi < m_jpsiMassLower || mass_Jpsi > m_jpsiMassUpper) {
531  ATH_MSG_DEBUG(" Original Jpsi candidate rejected by the mass cut: mass = "
532  << mass_Jpsi << " != (" << m_jpsiMassLower << ", " << m_jpsiMassUpper << ")" );
533  continue;
534  }
535 
536  // Check J/psi pi+ candidate invariant mass and skip if need be
537  double mass_Jpsipi = m_V0Tools->invariantMass(*vxcItr, massesJpsipi);
538  ATH_MSG_DEBUG("Jpsipi mass " << mass_Jpsipi);
539  if (mass_Jpsipi < m_jpsipiMassLower || mass_Jpsipi > m_jpsipiMassUpper) {
540  ATH_MSG_DEBUG(" Original Jpsipi candidate rejected by the mass cut: mass = "
541  << mass_Jpsipi << " != (" << m_jpsipiMassLower << ", " << m_jpsipiMassUpper << ")" );
542  continue;
543  }
544 
545  selectedJpsipiCandidates.push_back(*vxcItr);
546  }
547  if(selectedJpsipiCandidates.size()<1) return StatusCode::SUCCESS;
548 
549  // Select the D0/D0b candidates before calling cascade fit
550  std::vector<const xAOD::Vertex*> selectedD0Candidates;
551  for(auto vxcItr=d0Container->cbegin(); vxcItr!=d0Container->cend(); ++vxcItr) {
552 
553  // Check the passed flag first
554  const xAOD::Vertex* vtx = *vxcItr;
555  SG::AuxElement::Accessor<Char_t> flagAcc1("passed_D0");
556  SG::AuxElement::Accessor<Char_t> flagAcc2("passed_D0b");
557  bool isD0(true);
558  bool isD0b(true);
559  if(flagAcc1.isAvailable(*vtx)){
560  if(!flagAcc1(*vtx)) isD0 = false;
561  }
562  if(flagAcc2.isAvailable(*vtx)){
563  if(!flagAcc2(*vtx)) isD0b = false;
564  }
565  if(!(isD0||isD0b)) continue;
566 
567  // Ensure the total charge is correct
568  if ((*vxcItr)->trackParticle(0)->charge() != 1 || (*vxcItr)->trackParticle(1)->charge() != -1) {
569  ATH_MSG_DEBUG(" Original D0/D0-bar candidate rejected by the charge requirement: "
570  << (*vxcItr)->trackParticle(0)->charge() << ", " << (*vxcItr)->trackParticle(1)->charge() );
571  continue;
572  }
573 
574  // Check D0/D0bar candidate invariant mass and skip if need be
575  double mass_D0 = m_V0Tools->invariantMass(*vxcItr,massesD0);
576  double mass_D0b = m_V0Tools->invariantMass(*vxcItr,massesD0b);
577  ATH_MSG_DEBUG("D0 mass " << mass_D0 << ", D0b mass "<<mass_D0b);
578  if ((mass_D0 < m_D0MassLower || mass_D0 > m_D0MassUpper) && (mass_D0b < m_D0MassLower || mass_D0b > m_D0MassUpper)) {
579  ATH_MSG_DEBUG(" Original D0 candidate rejected by the mass cut: mass = "
580  << mass_D0 << " != (" << m_D0MassLower << ", " << m_D0MassUpper << ") "
581  << mass_D0b << " != (" << m_D0MassLower << ", " << m_D0MassUpper << ") " );
582  continue;
583  }
584 
585  selectedD0Candidates.push_back(*vxcItr);
586  }
587  if(selectedD0Candidates.size()<1) return StatusCode::SUCCESS;
588 
589  // Select J/psi D*+ candidates
590  // Iterate over Jpsi+pi vertices
591  for(auto jpsipiItr=selectedJpsipiCandidates.cbegin(); jpsipiItr!=selectedJpsipiCandidates.cend(); ++jpsipiItr) {
592 
593  size_t jpsipiTrkNum = (*jpsipiItr)->nTrackParticles();
594  tracksJpsipi.clear();
595  tracksJpsi.clear();
596  for( unsigned int it=0; it<jpsipiTrkNum; it++) tracksJpsipi.push_back((*jpsipiItr)->trackParticle(it));
597  for( unsigned int it=0; it<jpsipiTrkNum-1; it++) tracksJpsi.push_back((*jpsipiItr)->trackParticle(it));
598 
599  if (tracksJpsipi.size() != 3 || massesJpsipi.size() != 3 ) {
600  ATH_MSG_INFO("problems with Jpsi+pi input");
601  }
602 
603  bool tagD0(true);
604  if(abs(m_Dx_pid)==421 && (*jpsipiItr)->trackParticle(2)->charge()==-1) tagD0 = false;
605 
606  TLorentzVector p4_pi1; // Momentum of soft pion
607  p4_pi1.SetPtEtaPhiM((*jpsipiItr)->trackParticle(2)->pt(),
608  (*jpsipiItr)->trackParticle(2)->eta(),
609  (*jpsipiItr)->trackParticle(2)->phi(), m_vtx0Daug3MassHypo);
610 
611  // Iterate over D0/D0bar vertices
612  for(auto d0Itr=selectedD0Candidates.cbegin(); d0Itr!=selectedD0Candidates.cend(); ++d0Itr) {
613 
614  // Check identical tracks in input
615  if(std::find(tracksJpsipi.cbegin(), tracksJpsipi.cend(), (*d0Itr)->trackParticle(0)) != tracksJpsipi.cend()) continue;
616  if(std::find(tracksJpsipi.cbegin(), tracksJpsipi.cend(), (*d0Itr)->trackParticle(1)) != tracksJpsipi.cend()) continue;
617 
618 
619  TLorentzVector p4_ka, p4_pi2;
620  if(tagD0){ // for D*+
621  p4_pi2.SetPtEtaPhiM((*d0Itr)->trackParticle(0)->pt(),
622  (*d0Itr)->trackParticle(0)->eta(),
623  (*d0Itr)->trackParticle(0)->phi(), m_vtx1Daug1MassHypo);
624  p4_ka.SetPtEtaPhiM( (*d0Itr)->trackParticle(1)->pt(),
625  (*d0Itr)->trackParticle(1)->eta(),
626  (*d0Itr)->trackParticle(1)->phi(), m_vtx1Daug2MassHypo);
627  }else{ // change the order in the case of D*-
628  p4_pi2.SetPtEtaPhiM((*d0Itr)->trackParticle(1)->pt(),
629  (*d0Itr)->trackParticle(1)->eta(),
630  (*d0Itr)->trackParticle(1)->phi(), m_vtx1Daug1MassHypo);
631  p4_ka.SetPtEtaPhiM( (*d0Itr)->trackParticle(0)->pt(),
632  (*d0Itr)->trackParticle(0)->eta(),
633  (*d0Itr)->trackParticle(0)->phi(), m_vtx1Daug2MassHypo);
634  }
635  // Check D*+/- candidate invariant mass and skip if need be
636  double mass_Dst= (p4_pi1 + p4_ka + p4_pi2).M();
637  ATH_MSG_DEBUG("D*+/- mass " << mass_Dst);
638  if (mass_Dst < m_DstMassLower || mass_Dst > m_DstMassUpper) {
639  ATH_MSG_DEBUG(" Original D*+/- candidate rejected by the mass cut: mass = "
640  << mass_Dst << " != (" << m_DstMassLower << ", " << m_DstMassUpper << ")" );
641  continue;
642  }
643 
644  size_t d0TrkNum = (*d0Itr)->nTrackParticles();
645  tracksD0.clear();
646  for( unsigned int it=0; it<d0TrkNum; it++) tracksD0.push_back((*d0Itr)->trackParticle(it));
647  if (tracksD0.size() != 2 || massesD0.size() != 2 ) {
648  ATH_MSG_INFO("problems with D0 input");
649  }
650 
651  ATH_MSG_DEBUG("using tracks" << tracksJpsipi[0] << ", " << tracksJpsipi[1] << ", " << tracksJpsipi[2] << ", " << tracksD0[0] << ", " << tracksD0[1]);
652  ATH_MSG_DEBUG("Charge of Jpsi+pi tracks: "<<(*jpsipiItr)->trackParticle(0)->charge()<<", "<<(*jpsipiItr)->trackParticle(1)->charge()<<", "<<(*jpsipiItr)->trackParticle(2)->charge());
653  ATH_MSG_DEBUG("Charge of D0 tracks: "<<(*d0Itr)->trackParticle(0)->charge()<<", "<<(*d0Itr)->trackParticle(1)->charge());
654 
655  tracksBc.clear();
656  for( unsigned int it=0; it<jpsipiTrkNum; it++) tracksBc.push_back((*jpsipiItr)->trackParticle(it));
657  for( unsigned int it=0; it<d0TrkNum; it++) tracksBc.push_back((*d0Itr)->trackParticle(it));
658 
659 
660  // Apply the user's settings to the fitter
661  // Reset
662  std::unique_ptr<Trk::IVKalState> state (m_iVertexFitter->makeState());
663  // Robustness
664  int robustness = 0;
665  m_iVertexFitter->setRobustness(robustness, *state);
666  // Build up the topology
667  // Vertex list
668  std::vector<Trk::VertexID> vrtList;
669  // D0 vertex
670  Trk::VertexID vID;
671  if (m_constrD0) {
672  if(tagD0) vID = m_iVertexFitter->startVertex(tracksD0,massesD0,*state,mass_d0);
673  else vID = m_iVertexFitter->startVertex(tracksD0,massesD0b,*state,mass_d0);
674  } else {
675  if(tagD0) vID = m_iVertexFitter->startVertex(tracksD0,massesD0,*state);
676  else vID = m_iVertexFitter->startVertex(tracksD0,massesD0b,*state);
677  }
678  vrtList.push_back(vID);
679  // B vertex including Jpsi+pi
680  Trk::VertexID vID2 = m_iVertexFitter->nextVertex(tracksJpsipi,massesJpsipi,vrtList,*state);
681  if (m_constrJpsi) {
682  std::vector<Trk::VertexID> cnstV;
683  cnstV.clear();
684  if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJpsi,cnstV,*state,m_mass_jpsi).isSuccess() ) {
685  ATH_MSG_WARNING("addMassConstraint failed");
686  //return StatusCode::FAILURE;
687  }
688  }
689  // Do the work
690  std::unique_ptr<Trk::VxCascadeInfo> result(m_iVertexFitter->fitCascade(*state));
691 
692  if (result != nullptr) {
693 
694  // reset links to original tracks
695  BPhysPVCascadeTools::PrepareVertexLinks(result.get(), trackContainer);
696  ATH_MSG_DEBUG("storing tracks " << ((result->vertices())[0])->trackParticle(0) << ", "
697  << ((result->vertices())[0])->trackParticle(1) << ", "
698  << ((result->vertices())[1])->trackParticle(0) << ", "
699  << ((result->vertices())[1])->trackParticle(1) << ", "
700  << ((result->vertices())[1])->trackParticle(2));
701  // necessary to prevent memory leak
702  result->setSVOwnership(true);
703 
704  // Chi2/DOF cut
705  double bChi2DOF = result->fitChi2()/result->nDoF();
706  ATH_MSG_DEBUG("Candidate chi2/DOF is " << bChi2DOF);
707  bool chi2CutPassed = (m_chi2cut <= 0.0 || bChi2DOF < m_chi2cut);
708 
709  const std::vector< std::vector<TLorentzVector> > &moms = result->getParticleMoms();
710  double mass = m_CascadeTools->invariantMass(moms[1]);
711  if(chi2CutPassed) {
712  if (mass >= m_MassLower && mass <= m_MassUpper) {
713  cascadeinfoContainer->push_back(result.release());
714  } else {
715  ATH_MSG_DEBUG("Candidate rejected by the mass cut: mass = "
716  << mass << " != (" << m_MassLower << ", " << m_MassUpper << ")" );
717  }
718  }
719  }
720 
721  } //Iterate over D0 vertices
722 
723  } //Iterate over Jpsi+pi vertices
724 
725  ATH_MSG_DEBUG("cascadeinfoContainer size " << cascadeinfoContainer->size());
726 
727  return StatusCode::SUCCESS;
728  }
729 
730 }
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
DerivationFramework::JpsiPlusDpstCascade::JpsiPlusDpstCascade
JpsiPlusDpstCascade(const std::string &t, const std::string &n, const IInterface *p)
Definition: JpsiPlusDpstCascade.cxx:399
DerivationFramework::JpsiPlusDpstCascade::m_DstMassUpper
double m_DstMassUpper
Definition: JpsiPlusDpstCascade.h:60
xAOD::BPhysHypoHelper::setMass
bool setMass(const float val)
Set given invariant mass and its error.
Definition: BPhysHypoHelper.cxx:49
DerivationFramework::JpsiPlusDpstCascade::m_vtx0MassHypo
double m_vtx0MassHypo
Definition: JpsiPlusDpstCascade.h:63
xAOD::TrackParticle_v1::pt
virtual double pt() const override final
The transverse momentum ( ) of the particle.
Definition: TrackParticle_v1.cxx:73
DerivationFramework::JpsiPlusDpstCascade::m_PV_max
int m_PV_max
Definition: JpsiPlusDpstCascade.h:89
xAOD::Vertex_v1::x
float x() const
Returns the x position.
V0Tools.h
DerivationFramework::JpsiPlusDpstCascade::m_hypoName
std::string m_hypoName
name of the mass hypothesis.
Definition: JpsiPlusDpstCascade.h:86
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
DerivationFramework::JpsiPlusDpstCascade::m_MassUpper
double m_MassUpper
Definition: JpsiPlusDpstCascade.h:62
DerivationFramework::JpsiPlusDpstCascade::initialize
virtual StatusCode initialize() override
Definition: JpsiPlusDpstCascade.cxx:29
get_generator_info.result
result
Definition: get_generator_info.py:21
xAOD::VertexAuxContainer_v1
Temporary container used until we have I/O for AuxStoreInternal.
Definition: VertexAuxContainer_v1.h:32
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
Trk::VxCascadeInfo
Definition: VxCascadeInfo.h:75
DerivationFramework::JpsiPlusDpstCascade::m_vtx0Daug3MassHypo
double m_vtx0Daug3MassHypo
Definition: JpsiPlusDpstCascade.h:67
DerivationFramework::JpsiPlusDpstCascade::m_iVertexFitter
ToolHandle< Trk::TrkVKalVrtFitter > m_iVertexFitter
Definition: JpsiPlusDpstCascade.h:79
Trk::VertexID
int VertexID
Definition: IVertexCascadeFitter.h:23
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
JpsiPlusDpstCascade.h
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
SG::Accessor
Helper class to provide type-safe access to aux data.
Definition: Control/AthContainers/AthContainers/Accessor.h:66
xAOD::TrackParticle_v1::charge
float charge() const
Returns the charge.
Definition: TrackParticle_v1.cxx:150
DerivationFramework::JpsiPlusDpstCascade::m_DstMassLower
double m_DstMassLower
Definition: JpsiPlusDpstCascade.h:59
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::JpsiPlusDpstCascade::m_chi2cut
double m_chi2cut
Definition: JpsiPlusDpstCascade.h:76
skel.it
it
Definition: skel.GENtoEVGEN.py:423
test_pyathena.pt
pt
Definition: test_pyathena.py:11
DerivationFramework::JpsiPlusDpstCascade::m_V0Tools
ToolHandle< Trk::V0Tools > m_V0Tools
Definition: JpsiPlusDpstCascade.h:81
DerivationFramework::JpsiPlusDpstCascade::m_pvRefitter
ToolHandle< Analysis::PrimaryVertexRefitter > m_pvRefitter
Definition: JpsiPlusDpstCascade.h:80
DerivationFramework::JpsiPlusDpstCascade::m_PV_minNTracks
size_t m_PV_minNTracks
Definition: JpsiPlusDpstCascade.h:91
DataVector::cend
const_iterator cend() const noexcept
Return a const_iterator pointing past the end of the collection.
LArG4FSStartPointFilter.evt
evt
Definition: LArG4FSStartPointFilter.py:42
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
xAOD::VertexContainer
VertexContainer_v1 VertexContainer
Definition of the current "Vertex container version".
Definition: VertexContainer.h:14
x
#define x
dqt_zlumi_pandas.mass
mass
Definition: dqt_zlumi_pandas.py:170
DerivationFramework::JpsiPlusDpstCascade::m_DoVertexType
int m_DoVertexType
Definition: JpsiPlusDpstCascade.h:90
TrkVKalVrtFitter.h
DerivationFramework::JpsiPlusDpstCascade::m_VxPrimaryCandidateName
std::string m_VxPrimaryCandidateName
Name of primary vertex container.
Definition: JpsiPlusDpstCascade.h:51
DerivationFramework::JpsiPlusDpstCascade::m_vertexD0ContainerKey
std::string m_vertexD0ContainerKey
Definition: JpsiPlusDpstCascade.h:48
runBeamSpotCalibration.helper
helper
Definition: runBeamSpotCalibration.py:112
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::JpsiPlusDpstCascade::~JpsiPlusDpstCascade
~JpsiPlusDpstCascade()
Definition: JpsiPlusDpstCascade.cxx:469
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:58
xAOD::VertexAuxContainer
VertexAuxContainer_v1 VertexAuxContainer
Definition of the current jet auxiliary container.
Definition: VertexAuxContainer.h:19
lumiFormat.i
int i
Definition: lumiFormat.py:92
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::JpsiPlusDpstCascade::m_jpsiMassUpper
double m_jpsiMassUpper
Definition: JpsiPlusDpstCascade.h:54
DerivationFramework::BPhysPVCascadeTools::PrepareVertexLinks
static void PrepareVertexLinks(Trk::VxCascadeInfo *result, const xAOD::TrackParticleContainer *importedTrackCollection)
Definition: BPhysPVCascadeTools.cxx:204
MUON
xAOD::Muon MUON
D3PD INCLUDES.
Definition: TileCellFillerTool.h:37
DerivationFramework::JpsiPlusDpstCascade::m_jpsipiMassUpper
double m_jpsipiMassUpper
Definition: JpsiPlusDpstCascade.h:56
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
DerivationFramework::JpsiPlusDpstCascade::m_Dx_pid
int m_Dx_pid
Definition: JpsiPlusDpstCascade.h:73
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::JpsiPlusDpstCascade::m_vtx0Daug2MassHypo
double m_vtx0Daug2MassHypo
Definition: JpsiPlusDpstCascade.h:66
DerivationFramework::JpsiPlusDpstCascade::m_vtx1Daug1MassHypo
double m_vtx1Daug1MassHypo
Definition: JpsiPlusDpstCascade.h:68
DerivationFramework::VertexLink
ElementLink< xAOD::VertexContainer > VertexLink
Definition: Cascade3Plus1.cxx:23
DataVector
Derived DataVector<T>.
Definition: DataVector.h:581
BPhysPVCascadeTools.h
DerivationFramework::JpsiPlusDpstCascade::m_eventInfo_key
SG::ReadHandleKey< xAOD::EventInfo > m_eventInfo_key
Definition: JpsiPlusDpstCascade.h:78
DerivationFramework::JpsiPlusDpstCascade::m_refitPV
bool m_refitPV
Definition: JpsiPlusDpstCascade.h:84
DerivationFramework::JpsiPlusDpstCascade::m_cascadeOutputsKeys
std::vector< std::string > m_cascadeOutputsKeys
Definition: JpsiPlusDpstCascade.h:49
DerivationFramework::BPhysPVCascadeTools
Definition: BPhysPVCascadeTools.h:34
DerivationFramework::JpsiPlusDpstCascade::performSearch
StatusCode performSearch(std::vector< Trk::VxCascadeInfo * > *cascadeinfoContainer) const
Definition: JpsiPlusDpstCascade.cxx:471
CascadeTools.h
DerivationFramework::JpsiPlusDpstCascade::m_constrJpsi
bool m_constrJpsi
Definition: JpsiPlusDpstCascade.h:75
IVertexFitter.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
DerivationFramework::JpsiPlusDpstCascade::m_jpsiMassLower
double m_jpsiMassLower
Definition: JpsiPlusDpstCascade.h:53
VertexContainer.h
xAOD::Vertex_v1
Class describing a Vertex.
Definition: Vertex_v1.h:42
DerivationFramework::JpsiPlusDpstCascade::m_mass_jpsi
double m_mass_jpsi
Definition: JpsiPlusDpstCascade.h:72
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
xAOD::Vertex_v1::y
float y() const
Returns the y position.
DerivationFramework::JpsiPlusDpstCascade::m_vtx1Daug2MassHypo
double m_vtx1Daug2MassHypo
Definition: JpsiPlusDpstCascade.h:69
declareProperty
#define declareProperty(n, p, h)
Definition: BaseFakeBkgTool.cxx:15
SG::ConstAccessor< T, AuxAllocator_t< T > >::isAvailable
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
DerivationFramework::JpsiPlusDpstCascade::m_refPVContainerName
std::string m_refPVContainerName
Definition: JpsiPlusDpstCascade.h:85
DerivationFramework::JpsiPlusDpstCascade::m_vtx1MassHypo
double m_vtx1MassHypo
Definition: JpsiPlusDpstCascade.h:64
DataVector::cbegin
const_iterator cbegin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
DerivationFramework::JpsiPlusDpstCascade::m_constrD0
bool m_constrD0
Definition: JpsiPlusDpstCascade.h:74
DerivationFramework::JpsiPlusDpstCascade::m_MassLower
double m_MassLower
Definition: JpsiPlusDpstCascade.h:61
DerivationFramework::JpsiPlusDpstCascade::addBranches
virtual StatusCode addBranches() const override
Pass the thinning service
Definition: JpsiPlusDpstCascade.cxx:62
AthAlgTool
Definition: AthAlgTool.h:26
DataVector::size
size_type size() const noexcept
Returns the number of elements in the collection.
DerivationFramework::JpsiPlusDpstCascade::m_vertexContainerKey
std::string m_vertexContainerKey
Definition: JpsiPlusDpstCascade.h:47
DerivationFramework::JpsiPlusDpstCascade::m_D0MassLower
double m_D0MassLower
Definition: JpsiPlusDpstCascade.h:57
HepMCHelpers.h
xAOD::BPhysHypoHelper::setMassErr
bool setMassErr(const float val)
invariant mass error
Definition: BPhysHypoHelper.cxx:54
DerivationFramework::JpsiPlusDpstCascade::m_D0MassUpper
double m_D0MassUpper
Definition: JpsiPlusDpstCascade.h:58
DerivationFramework::JpsiPlusDpstCascade::m_jpsipiMassLower
double m_jpsipiMassLower
Definition: JpsiPlusDpstCascade.h:55
VertexAuxContainer.h
xAOD::TrackParticle_v1::phi
virtual double phi() const override final
The azimuthal angle ( ) of the particle (has range to .)
DerivationFramework::JpsiPlusDpstCascade::m_CascadeTools
ToolHandle< DerivationFramework::CascadeTools > m_CascadeTools
Definition: JpsiPlusDpstCascade.h:82
DerivationFramework::JpsiPlusDpstCascade::m_vtx0Daug1MassHypo
double m_vtx0Daug1MassHypo
Definition: JpsiPlusDpstCascade.h:65
Generated on Tue Jul 2 2024 21:19:07 for ATLAS Offline Software by   doxygen 1.8.18