ATLAS Offline Software
Loading...
Searching...
No Matches
JpsiPlusDsCascade.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3*/
5// JpsiPlusDsCascade.cxx, (c) ATLAS Detector software
7#include "JpsiPlusDsCascade.h"
12#include "BPhysPVCascadeTools.h"
15#include <algorithm>
17#include "HepPDT/ParticleDataTable.hh"
18
21
22namespace DerivationFramework {
23 typedef ElementLink<xAOD::VertexContainer> VertexLink;
24 typedef std::vector<VertexLink> VertexLinkVector;
25 typedef std::vector<const xAOD::TrackParticle*> TrackBag;
26
28
29 // retrieving vertex Fitter
30 ATH_CHECK( m_iVertexFitter.retrieve());
31
32 // retrieving the V0 tools
33 ATH_CHECK( m_V0Tools.retrieve());
34
35 // retrieving the Cascade tools
36 ATH_CHECK( m_CascadeTools.retrieve());
37
38 // Get the beam spot service
39 ATH_CHECK( m_eventInfo_key.initialize() );
40
41 ATH_CHECK( m_partPropSvc.retrieve() );
42 auto pdt = m_partPropSvc->PDT();
43
44 // retrieve particle masses
46 if(m_vtx0MassHypo < 0.)
48 if(m_vtx1MassHypo < 0.) {
51 }
52
55 if(m_vtx1Daug1MassHypo < 0.) {
58 }
61
62 return StatusCode::SUCCESS;
63 }
64
65
66 StatusCode JpsiPlusDsCascade::addBranches(const EventContext& ctx) const
67 {
68 std::vector<Trk::VxCascadeInfo*> cascadeinfoContainer;
69 constexpr int topoN = 2;
70 std::array<xAOD::VertexContainer*, topoN> Vtxwritehandles;
71 std::array<xAOD::VertexAuxContainer*, topoN> Vtxwritehandlesaux;
72 if(m_cascadeOutputsKeys.size() !=topoN) { ATH_MSG_FATAL("Incorrect number of VtxContainers"); return StatusCode::FAILURE; }
73
74 for(int i =0; i<topoN;i++){
75 Vtxwritehandles[i] = new xAOD::VertexContainer();
76 Vtxwritehandlesaux[i] = new xAOD::VertexAuxContainer();
77 Vtxwritehandles[i]->setStore(Vtxwritehandlesaux[i]);
78 ATH_CHECK(evtStore()->record(Vtxwritehandles[i] , m_cascadeOutputsKeys[i] )); // FIXME Use Handles
79 ATH_CHECK(evtStore()->record(Vtxwritehandlesaux[i], m_cascadeOutputsKeys[i] + "Aux.")); // FIXME Use Handles
80 }
81
82 //----------------------------------------------------
83 // retrieve primary vertices
84 //----------------------------------------------------
85 const xAOD::Vertex * primaryVertex{};
86 const xAOD::VertexContainer *pvContainer{};
87 ATH_CHECK(evtStore()->retrieve(pvContainer, m_VxPrimaryCandidateName)); // FIXME Use Handles
88 ATH_MSG_DEBUG("Found " << m_VxPrimaryCandidateName << " in StoreGate!");
89
90 if (pvContainer->size()==0){
91 ATH_MSG_WARNING("You have no primary vertices: " << pvContainer->size());
92 return StatusCode::RECOVERABLE;
93 } else {
94 primaryVertex = (*pvContainer)[0];
95 }
96
97 //----------------------------------------------------
98 // Try to retrieve refitted primary vertices
99 //----------------------------------------------------
100 xAOD::VertexContainer* refPvContainer{};
101 xAOD::VertexAuxContainer* refPvAuxContainer{};
102 if (m_refitPV) {
104 // refitted PV container exists. Get it from the store gate
105 ATH_CHECK(evtStore()->retrieve(refPvContainer , m_refPVContainerName )); // FIXME Use Handles
106 ATH_CHECK(evtStore()->retrieve(refPvAuxContainer, m_refPVContainerName + "Aux.")); // FIXME Use Handles
107 } else {
108 // refitted PV container does not exist. Create a new one.
109 refPvContainer = new xAOD::VertexContainer;
110 refPvAuxContainer = new xAOD::VertexAuxContainer;
111 refPvContainer->setStore(refPvAuxContainer);
112 ATH_CHECK(evtStore()->record(refPvContainer , m_refPVContainerName)); // FIXME Use Handles
113 ATH_CHECK(evtStore()->record(refPvAuxContainer, m_refPVContainerName+"Aux.")); // FIXME Use Handles
114 }
115 }
116
117 ATH_CHECK(performSearch(&cascadeinfoContainer, ctx));
118
120 if(not evt.isValid()) ATH_MSG_ERROR("Cannot Retrieve " << m_eventInfo_key.key() );
121 BPhysPVCascadeTools helper(&(*m_CascadeTools), evt.cptr());
122 helper.SetMinNTracksInPV(m_PV_minNTracks);
123
124 // Decorators for the main vertex: chi2, ndf, pt and pt error, plus the V0 vertex variables
125 SG::AuxElement::Decorator<VertexLinkVector> CascadeLinksDecor("CascadeVertexLinks");
126 SG::AuxElement::Decorator<VertexLinkVector> JpsiLinksDecor("JpsiVertexLinks");
127 SG::AuxElement::Decorator<VertexLinkVector> DxLinksDecor("DxVertexLinks");
128 SG::AuxElement::Decorator<float> chi2_decor("ChiSquared");
129 SG::AuxElement::Decorator<float> ndof_decor("NumberDoF");
131 SG::AuxElement::Decorator<float> PtErr_decor("PtErr");
132 SG::AuxElement::Decorator<float> Mass_svdecor("Dx_mass");
133 SG::AuxElement::Decorator<float> MassErr_svdecor("Dx_massErr");
134 SG::AuxElement::Decorator<float> Pt_svdecor("Dx_Pt");
135 SG::AuxElement::Decorator<float> PtErr_svdecor("Dx_PtErr");
136 SG::AuxElement::Decorator<float> Lxy_svdecor("Dx_Lxy");
137 SG::AuxElement::Decorator<float> LxyErr_svdecor("Dx_LxyErr");
138 SG::AuxElement::Decorator<float> Tau_svdecor("Dx_Tau");
139 SG::AuxElement::Decorator<float> TauErr_svdecor("Dx_TauErr");
140
141 SG::AuxElement::Decorator<float> MassMumu_decor("Mumu_mass");
142 SG::AuxElement::Decorator<float> MassKX_svdecor("KX_mass");
143 SG::AuxElement::Decorator<float> MassKXpi_svdecor("KXpi_mass");
144
145 ATH_MSG_DEBUG("cascadeinfoContainer size " << cascadeinfoContainer.size());
146
147 // Get Jpsi container and identify the input Jpsi
148 const xAOD::VertexContainer *jpsiContainer{};
149 ATH_CHECK(evtStore()->retrieve(jpsiContainer , m_vertexContainerKey )); // FIXME Use Handles
150 const xAOD::VertexContainer *dxContainer{};
151 ATH_CHECK(evtStore()->retrieve(dxContainer , m_vertexDxContainerKey )); // FIXME Use Handles
152
153 for (Trk::VxCascadeInfo* x : cascadeinfoContainer) {
154 if(x==nullptr) {
155 ATH_MSG_ERROR("cascadeinfoContainer is null");
156 //x is dereferenced if we pass this
157 return StatusCode::FAILURE;
158 }
159
160 // the cascade fitter returns:
161 // std::vector<xAOD::Vertex*>, each xAOD::Vertex contains the refitted track parameters (perigee at the vertex position)
162 // vertices[iv] the links to the original TPs and a covariance of size 3+5*NTRK; the chi2 of the total fit
163 // is split between the cascade vertices as per track contribution
164 // std::vector< std::vector<TLorentzVector> >, each std::vector<TLorentzVector> contains the refitted momenta (TLorentzVector)
165 // momenta[iv][...] of all tracks in the corresponding vertex, including any pseudotracks (from cascade vertices)
166 // originating in this vertex; the masses are as assigned in the cascade fit
167 // std::vector<Amg::MatrixX>, the corresponding covariance matrices in momentum space
168 // covariance[iv]
169 // int nDoF, double Chi2
170 //
171 // the invariant mass, pt, lifetime etc. errors should be calculated using the covariance matrices in momentum space as these
172 // take into account the full track-track and track-vertex correlations
173 //
174 // 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.
175 // The covariance terms between the two vertices are not stored. In momentum space momenta[0] contains the 2 V0 tracks,
176 // their momenta add up to the momentum of the 3rd track in momenta[1], the first two being the Jpsi tracks
177
178 const std::vector<xAOD::Vertex*> &cascadeVertices = x->vertices();
179 if(cascadeVertices.size()!=topoN)
180 ATH_MSG_ERROR("Incorrect number of vertices");
181 if(cascadeVertices[0] == nullptr || cascadeVertices[1] == nullptr) ATH_MSG_ERROR("Error null vertex");
182 // Keep vertices (bear in mind that they come in reverse order!)
183 for(int i =0;i<topoN;i++) Vtxwritehandles[i]->push_back(cascadeVertices[i]);
184
185 x->setSVOwnership(false); // Prevent Container from deleting vertices
186 const auto mainVertex = cascadeVertices[1]; // this is the B_c+/- vertex
187 const std::vector< std::vector<TLorentzVector> > &moms = x->getParticleMoms();
188
189 // Set links to cascade vertices
190 std::vector<const xAOD::Vertex*> verticestoLink;
191 verticestoLink.push_back(cascadeVertices[0]);
192 if(Vtxwritehandles[1] == nullptr) ATH_MSG_ERROR("Vtxwritehandles[1] is null");
193 if(!BPhysPVCascadeTools::LinkVertices(CascadeLinksDecor, verticestoLink, Vtxwritehandles[0], cascadeVertices[1]))
194 ATH_MSG_ERROR("Error decorating with cascade vertices");
195
196 // Identify the input Jpsi
197 const xAOD::Vertex* jpsiVertex = BPhysPVCascadeTools::FindVertex<2>(jpsiContainer, cascadeVertices[1]);
198 ATH_MSG_DEBUG("1 pt Jpsi tracks " << cascadeVertices[1]->trackParticle(0)->pt() << ", " << cascadeVertices[1]->trackParticle(1)->pt());
199 if (jpsiVertex) ATH_MSG_DEBUG("2 pt Jpsi tracks " << jpsiVertex->trackParticle(0)->pt() << ", " << jpsiVertex->trackParticle(1)->pt());
200
201 // Identify the input D_(s)+
202 const xAOD::Vertex* dxVertex = BPhysPVCascadeTools::FindVertex<3>(dxContainer, cascadeVertices[0]);;
203 ATH_MSG_DEBUG("1 pt D_(s)+ tracks " << cascadeVertices[0]->trackParticle(0)->pt() << ", " << cascadeVertices[0]->trackParticle(1)->pt() << ", " << cascadeVertices[0]->trackParticle(2)->pt());
204 if (dxVertex) ATH_MSG_DEBUG("2 pt D_(s)+ tracks " << dxVertex->trackParticle(0)->pt() << ", " << dxVertex->trackParticle(1)->pt() << ", " << dxVertex->trackParticle(2)->pt());
205
206 // Set links to input vertices
207 std::vector<const xAOD::Vertex*> jpsiVerticestoLink;
208 if (jpsiVertex) jpsiVerticestoLink.push_back(jpsiVertex);
209 else ATH_MSG_WARNING("Could not find linking Jpsi");
210 if(!BPhysPVCascadeTools::LinkVertices(JpsiLinksDecor, jpsiVerticestoLink, jpsiContainer, cascadeVertices[1]))
211 ATH_MSG_ERROR("Error decorating with Jpsi vertices");
212
213 std::vector<const xAOD::Vertex*> dxVerticestoLink;
214 if (dxVertex) dxVerticestoLink.push_back(dxVertex);
215 else ATH_MSG_WARNING("Could not find linking D_(s)+");
216 if(!BPhysPVCascadeTools::LinkVertices(DxLinksDecor, dxVerticestoLink, dxContainer, cascadeVertices[1]))
217 ATH_MSG_ERROR("Error decorating with D_(s)+ vertices");
218
219 bool tagDp(true);
220 if (dxVertex) {
221 if(abs(m_Dx_pid)==411 && (dxVertex->trackParticle(2)->charge()==-1)) tagDp = false;
222 }
223
224 double mass_b = m_vtx0MassHypo;
225 double mass_d = m_vtx1MassHypo;
226 std::vector<double> massesJpsi;
227 massesJpsi.push_back(m_vtx0Daug1MassHypo);
228 massesJpsi.push_back(m_vtx0Daug2MassHypo);
229 std::vector<double> massesDx;
230 if(tagDp){
231 massesDx.push_back(m_vtx1Daug1MassHypo);
232 massesDx.push_back(m_vtx1Daug2MassHypo);
233 }else{ // Change the order for D-
234 massesDx.push_back(m_vtx1Daug2MassHypo);
235 massesDx.push_back(m_vtx1Daug1MassHypo);
236 }
237 massesDx.push_back(m_vtx1Daug3MassHypo);
238 std::vector<double> Masses;
239 Masses.push_back(m_vtx0Daug1MassHypo);
240 Masses.push_back(m_vtx0Daug2MassHypo);
241 Masses.push_back(m_vtx1MassHypo);
242
243 // loop over candidates -- Don't apply PV_minNTracks requirement here
244 // because it may result in exclusion of the high-pt PV.
245 // get good PVs
246
247
248 xAOD::BPhysHypoHelper vtx(m_hypoName, mainVertex);
249
250 // Get refitted track momenta from all vertices, charged tracks only
252
253 // Decorate main vertex
254 //
255 // 1.a) mass, mass error
256 BPHYS_CHECK( vtx.setMass(m_CascadeTools->invariantMass(moms[1])) );
257 BPHYS_CHECK( vtx.setMassErr(m_CascadeTools->invariantMassError(moms[1],x->getCovariance()[1])) );
258 // 1.b) pt and pT error (the default pt of mainVertex is != the pt of the full cascade fit!)
259 Pt_decor(*mainVertex) = m_CascadeTools->pT(moms[1]);
260 PtErr_decor(*mainVertex) = m_CascadeTools->pTError(moms[1],x->getCovariance()[1]);
261 // 1.c) chi2 and ndof (the default chi2 of mainVertex is != the chi2 of the full cascade fit!)
262 chi2_decor(*mainVertex) = x->fitChi2();
263 ndof_decor(*mainVertex) = x->nDoF();
264
265 float massMumu = 0.;
266 if (jpsiVertex) {
267 TLorentzVector p4_mu1, p4_mu2;
268 p4_mu1.SetPtEtaPhiM(jpsiVertex->trackParticle(0)->pt(),
269 jpsiVertex->trackParticle(0)->eta(),
270 jpsiVertex->trackParticle(0)->phi(), m_vtx0Daug1MassHypo);
271 p4_mu2.SetPtEtaPhiM(jpsiVertex->trackParticle(1)->pt(),
272 jpsiVertex->trackParticle(1)->eta(),
273 jpsiVertex->trackParticle(1)->phi(), m_vtx0Daug2MassHypo);
274 massMumu = (p4_mu1 + p4_mu2).M();
275 }
276 MassMumu_decor(*mainVertex) = massMumu;
277
278 float massKX = 0., massKXpi = 0.;
279 if (dxVertex) {
280 TLorentzVector p4_h1, p4_h2, p4_h3;
281 if(tagDp){
282 p4_h1.SetPtEtaPhiM(dxVertex->trackParticle(0)->pt(),
283 dxVertex->trackParticle(0)->eta(),
284 dxVertex->trackParticle(0)->phi(), m_vtx1Daug1MassHypo);
285 p4_h2.SetPtEtaPhiM(dxVertex->trackParticle(1)->pt(),
286 dxVertex->trackParticle(1)->eta(),
287 dxVertex->trackParticle(1)->phi(), m_vtx1Daug2MassHypo);
288 }else{ // Change the order for D-
289 p4_h1.SetPtEtaPhiM(dxVertex->trackParticle(0)->pt(),
290 dxVertex->trackParticle(0)->eta(),
291 dxVertex->trackParticle(0)->phi(), m_vtx1Daug2MassHypo);
292 p4_h2.SetPtEtaPhiM(dxVertex->trackParticle(1)->pt(),
293 dxVertex->trackParticle(1)->eta(),
294 dxVertex->trackParticle(1)->phi(), m_vtx1Daug1MassHypo);
295 }
296 p4_h3.SetPtEtaPhiM(dxVertex->trackParticle(2)->pt(),
297 dxVertex->trackParticle(2)->eta(),
298 dxVertex->trackParticle(2)->phi(), m_vtx1Daug3MassHypo);
299 massKX = (p4_h1 + p4_h2).M();
300 massKXpi = (p4_h1 + p4_h2 + p4_h3).M();
301 }
302 MassKX_svdecor(*mainVertex) = massKX;
303 MassKXpi_svdecor(*mainVertex) = massKXpi;
304
305 ATH_CHECK(helper.FillCandwithRefittedVertices(m_refitPV, pvContainer,
306 refPvContainer, &(*m_pvRefitter), m_PV_max, m_DoVertexType, x, 1, mass_b, vtx));
307
308
309 // 4) decorate the main vertex with V0 vertex mass, pt, lifetime and lxy values (plus errors)
310 // V0 points to the main vertex, so lifetime and lxy are w.r.t the main vertex
311 Mass_svdecor(*mainVertex) = m_CascadeTools->invariantMass(moms[0]);
312 MassErr_svdecor(*mainVertex) = m_CascadeTools->invariantMassError(moms[0],x->getCovariance()[0]);
313 Pt_svdecor(*mainVertex) = m_CascadeTools->pT(moms[0]);
314 PtErr_svdecor(*mainVertex) = m_CascadeTools->pTError(moms[0],x->getCovariance()[0]);
315 Lxy_svdecor(*mainVertex) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
316 LxyErr_svdecor(*mainVertex) = m_CascadeTools->lxyError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
317 Tau_svdecor(*mainVertex) = m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1]);
318 TauErr_svdecor(*mainVertex) = m_CascadeTools->tauError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
319
320 // Some checks in DEBUG mode
321 ATH_MSG_DEBUG("chi2 " << x->fitChi2()
322 << " chi2_1 " << m_V0Tools->chisq(cascadeVertices[0])
323 << " chi2_2 " << m_V0Tools->chisq(cascadeVertices[1])
324 << " vprob " << m_CascadeTools->vertexProbability(x->nDoF(),x->fitChi2()));
325 ATH_MSG_DEBUG("ndf " << x->nDoF() << " ndf_1 " << m_V0Tools->ndof(cascadeVertices[0]) << " ndf_2 " << m_V0Tools->ndof(cascadeVertices[1]));
326 ATH_MSG_DEBUG("V0Tools mass_d " << m_V0Tools->invariantMass(cascadeVertices[0],massesDx)
327 << " error " << m_V0Tools->invariantMassError(cascadeVertices[0],massesDx)
328 << " mass_J " << m_V0Tools->invariantMass(cascadeVertices[1],massesJpsi)
329 << " error " << m_V0Tools->invariantMassError(cascadeVertices[1],massesJpsi));
330 // masses and errors, using track masses assigned in the fit
331 double Mass_B = m_CascadeTools->invariantMass(moms[1]);
332 double Mass_D = m_CascadeTools->invariantMass(moms[0]);
333 double Mass_B_err = m_CascadeTools->invariantMassError(moms[1],x->getCovariance()[1]);
334 double Mass_D_err = m_CascadeTools->invariantMassError(moms[0],x->getCovariance()[0]);
335 ATH_MSG_DEBUG("Mass_B " << Mass_B << " Mass_D " << Mass_D);
336 ATH_MSG_DEBUG("Mass_B_err " << Mass_B_err << " Mass_D_err " << Mass_D_err);
337 double mprob_B = m_CascadeTools->massProbability(mass_b,Mass_B,Mass_B_err);
338 double mprob_D = m_CascadeTools->massProbability(mass_d,Mass_D,Mass_D_err);
339 ATH_MSG_DEBUG("mprob_B " << mprob_B << " mprob_D " << mprob_D);
340 // masses and errors, assigning user defined track masses
341 ATH_MSG_DEBUG("Mass_b " << m_CascadeTools->invariantMass(moms[1],Masses)
342 << " Mass_d " << m_CascadeTools->invariantMass(moms[0],massesDx));
343 ATH_MSG_DEBUG("Mass_b_err " << m_CascadeTools->invariantMassError(moms[1],x->getCovariance()[1],Masses)
344 << " Mass_d_err " << m_CascadeTools->invariantMassError(moms[0],x->getCovariance()[0],massesDx));
345 ATH_MSG_DEBUG("pt_b " << m_CascadeTools->pT(moms[1])
346 << " pt_d " << m_CascadeTools->pT(moms[0])
347 << " pt_dp " << m_V0Tools->pT(cascadeVertices[0]));
348 ATH_MSG_DEBUG("ptErr_b " << m_CascadeTools->pTError(moms[1],x->getCovariance()[1])
349 << " ptErr_d " << m_CascadeTools->pTError(moms[0],x->getCovariance()[0])
350 << " ptErr_dp " << m_V0Tools->pTError(cascadeVertices[0]));
351 ATH_MSG_DEBUG("lxy_B " << m_V0Tools->lxy(cascadeVertices[1],primaryVertex) << " lxy_D " << m_V0Tools->lxy(cascadeVertices[0],cascadeVertices[1]));
352 ATH_MSG_DEBUG("lxy_b " << m_CascadeTools->lxy(moms[1],cascadeVertices[1],primaryVertex) << " lxy_d " << m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]));
353 ATH_MSG_DEBUG("lxyErr_b " << m_CascadeTools->lxyError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
354 << " lxyErr_d " << m_CascadeTools->lxyError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
355 << " lxyErr_dp " << m_V0Tools->lxyError(cascadeVertices[0],cascadeVertices[1]));
356 ATH_MSG_DEBUG("tau_B " << m_CascadeTools->tau(moms[1],cascadeVertices[1],primaryVertex,mass_b)
357 << " tau_dp " << m_V0Tools->tau(cascadeVertices[0],cascadeVertices[1],massesDx));
358 ATH_MSG_DEBUG("tau_b " << m_CascadeTools->tau(moms[1],cascadeVertices[1],primaryVertex)
359 << " tau_d " << m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1])
360 << " tau_D " << m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1],mass_d));
361 ATH_MSG_DEBUG("tauErr_b " << m_CascadeTools->tauError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
362 << " tauErr_d " << m_CascadeTools->tauError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
363 << " tauErr_dp " << m_V0Tools->tauError(cascadeVertices[0],cascadeVertices[1],massesDx));
364 ATH_MSG_DEBUG("TauErr_b " << m_CascadeTools->tauError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex,mass_b)
365 << " TauErr_d " << m_CascadeTools->tauError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1],mass_d)
366 << " TauErr_dp " << m_V0Tools->tauError(cascadeVertices[0],cascadeVertices[1],massesDx,mass_d));
367
368 ATH_MSG_DEBUG("CascadeTools main vert wrt PV " << " CascadeTools SV " << " V0Tools SV");
369 ATH_MSG_DEBUG("a0z " << m_CascadeTools->a0z(moms[1],cascadeVertices[1],primaryVertex)
370 << ", " << m_CascadeTools->a0z(moms[0],cascadeVertices[0],cascadeVertices[1])
371 << ", " << m_V0Tools->a0z(cascadeVertices[0],cascadeVertices[1]));
372 ATH_MSG_DEBUG("a0zErr " << m_CascadeTools->a0zError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
373 << ", " << m_CascadeTools->a0zError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
374 << ", " << m_V0Tools->a0zError(cascadeVertices[0],cascadeVertices[1]));
375 ATH_MSG_DEBUG("a0xy " << m_CascadeTools->a0xy(moms[1],cascadeVertices[1],primaryVertex)
376 << ", " << m_CascadeTools->a0xy(moms[0],cascadeVertices[0],cascadeVertices[1])
377 << ", " << m_V0Tools->a0xy(cascadeVertices[0],cascadeVertices[1]));
378 ATH_MSG_DEBUG("a0xyErr " << m_CascadeTools->a0xyError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
379 << ", " << m_CascadeTools->a0xyError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
380 << ", " << m_V0Tools->a0xyError(cascadeVertices[0],cascadeVertices[1]));
381 ATH_MSG_DEBUG("a0 " << m_CascadeTools->a0(moms[1],cascadeVertices[1],primaryVertex)
382 << ", " << m_CascadeTools->a0(moms[0],cascadeVertices[0],cascadeVertices[1])
383 << ", " << m_V0Tools->a0(cascadeVertices[0],cascadeVertices[1]));
384 ATH_MSG_DEBUG("a0Err " << m_CascadeTools->a0Error(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
385 << ", " << m_CascadeTools->a0Error(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
386 << ", " << m_V0Tools->a0Error(cascadeVertices[0],cascadeVertices[1]));
387 ATH_MSG_DEBUG("x0 " << m_V0Tools->vtx(cascadeVertices[0]).x() << " y0 " << m_V0Tools->vtx(cascadeVertices[0]).y() << " z0 " << m_V0Tools->vtx(cascadeVertices[0]).z());
388 ATH_MSG_DEBUG("x1 " << m_V0Tools->vtx(cascadeVertices[1]).x() << " y1 " << m_V0Tools->vtx(cascadeVertices[1]).y() << " z1 " << m_V0Tools->vtx(cascadeVertices[1]).z());
389 ATH_MSG_DEBUG("X0 " << primaryVertex->x() << " Y0 " << primaryVertex->y() << " Z0 " << primaryVertex->z());
390 ATH_MSG_DEBUG("rxy0 " << m_V0Tools->rxy(cascadeVertices[0]) << " rxyErr0 " << m_V0Tools->rxyError(cascadeVertices[0]));
391 ATH_MSG_DEBUG("rxy1 " << m_V0Tools->rxy(cascadeVertices[1]) << " rxyErr1 " << m_V0Tools->rxyError(cascadeVertices[1]));
392 ATH_MSG_DEBUG("Rxy0 wrt PV " << m_V0Tools->rxy(cascadeVertices[0],primaryVertex) << " RxyErr0 wrt PV " << m_V0Tools->rxyError(cascadeVertices[0],primaryVertex));
393 ATH_MSG_DEBUG("Rxy1 wrt PV " << m_V0Tools->rxy(cascadeVertices[1],primaryVertex) << " RxyErr1 wrt PV " << m_V0Tools->rxyError(cascadeVertices[1],primaryVertex));
394 ATH_MSG_DEBUG("number of covariance matrices " << (x->getCovariance()).size());
395 } // loop over cascadeinfoContainer
396
397 // Deleting cascadeinfo since this won't be stored.
398 // Vertices have been kept in m_cascadeOutputs and should be owned by their container
399 for (auto x : cascadeinfoContainer) delete x;
400
401 return StatusCode::SUCCESS;
402 }
403
404
405 JpsiPlusDsCascade::JpsiPlusDsCascade(const std::string& t, const std::string& n, const IInterface* p) : base_class(t,n,p),
408 m_cascadeOutputsKeys{ "JpsiPlusDsCascadeVtx1", "JpsiPlusDsCascadeVtx2" },
409 m_VxPrimaryCandidateName("PrimaryVertices"),
410 m_jpsiMassLower(0.0),
411 m_jpsiMassUpper(10000.0),
412 m_DxMassLower(0.0),
413 m_DxMassUpper(10000.0),
414 m_MassLower(0.0),
415 m_MassUpper(20000.0),
416 m_vtx0MassHypo(-1),
417 m_vtx1MassHypo(-1),
423 m_mass_jpsi(-1),
424 m_Dx_pid(431),
425 m_constrDx(true),
426 m_constrJpsi(true),
427 m_chi2cut(-1.0),
428 m_iVertexFitter("Trk::TrkVKalVrtFitter"),
429 m_pvRefitter("Analysis::PrimaryVertexRefitter"),
430 m_V0Tools("Trk::V0Tools"),
431 m_CascadeTools("DerivationFramework::CascadeTools")
432 {
433 declareProperty("JpsiVertices", m_vertexContainerKey);
434 declareProperty("DxVertices", m_vertexDxContainerKey);
435 declareProperty("VxPrimaryCandidateName", m_VxPrimaryCandidateName);
436 declareProperty("RefPVContainerName", m_refPVContainerName = "RefittedPrimaryVertices");
437 declareProperty("JpsiMassLowerCut", m_jpsiMassLower);
438 declareProperty("JpsiMassUpperCut", m_jpsiMassUpper);
439 declareProperty("DxMassLowerCut", m_DxMassLower);
440 declareProperty("DxMassUpperCut", m_DxMassUpper);
441 declareProperty("MassLowerCut", m_MassLower);
442 declareProperty("MassUpperCut", m_MassUpper);
443 declareProperty("HypothesisName", m_hypoName = "Bc");
444 declareProperty("Vtx0MassHypo", m_vtx0MassHypo);
445 declareProperty("Vtx1MassHypo", m_vtx1MassHypo);
446 declareProperty("Vtx0Daug1MassHypo", m_vtx0Daug1MassHypo);
447 declareProperty("Vtx0Daug2MassHypo", m_vtx0Daug2MassHypo);
448 declareProperty("Vtx1Daug1MassHypo", m_vtx1Daug1MassHypo);
449 declareProperty("Vtx1Daug2MassHypo", m_vtx1Daug2MassHypo);
450 declareProperty("Vtx1Daug3MassHypo", m_vtx1Daug3MassHypo);
451 declareProperty("JpsiMass", m_mass_jpsi);
452 declareProperty("DxHypothesis", m_Dx_pid);
453 declareProperty("ApplyDxMassConstraint", m_constrDx);
454 declareProperty("ApplyJpsiMassConstraint", m_constrJpsi);
455 declareProperty("Chi2Cut", m_chi2cut);
456 declareProperty("RefitPV", m_refitPV = true);
457 declareProperty("MaxnPV", m_PV_max = 999);
458 declareProperty("MinNTracksInPV", m_PV_minNTracks = 0);
459 declareProperty("DoVertexType", m_DoVertexType = 7);
460 declareProperty("TrkVertexFitterTool", m_iVertexFitter);
461 declareProperty("PVRefitter", m_pvRefitter);
462 declareProperty("V0Tools", m_V0Tools);
463 declareProperty("CascadeTools", m_CascadeTools);
464 declareProperty("CascadeVertexCollections", m_cascadeOutputsKeys);
465 }
466
468
469 StatusCode JpsiPlusDsCascade::performSearch(std::vector<Trk::VxCascadeInfo*> *cascadeinfoContainer, const EventContext&) const
470 {
471 ATH_MSG_DEBUG( "JpsiPlusDsCascade::performSearch" );
472 assert(cascadeinfoContainer!=nullptr);
473
474 // Get TrackParticle container (for setting links to the original tracks)
475 const xAOD::TrackParticleContainer *trackContainer{};
476 ATH_CHECK(evtStore()->retrieve(trackContainer , "InDetTrackParticles" )); // FIXME Use Handles
477
478 // Get Jpsi container
479 const xAOD::VertexContainer *jpsiContainer{};
480 ATH_CHECK(evtStore()->retrieve(jpsiContainer , m_vertexContainerKey )); // FIXME Use Handles
481
482 // Get V0 container
483 const xAOD::VertexContainer *dxContainer{};
484 ATH_CHECK(evtStore()->retrieve(dxContainer , m_vertexDxContainerKey )); // FIXME Use Handles
485
486
487
488 double mass_d = m_vtx1MassHypo;
489 std::vector<const xAOD::TrackParticle*> tracksJpsi;
490 std::vector<const xAOD::TrackParticle*> tracksDx;
491 std::vector<const xAOD::TrackParticle*> tracksBc;
492 std::vector<double> massesJpsi;
493 massesJpsi.push_back(m_vtx0Daug1MassHypo);
494 massesJpsi.push_back(m_vtx0Daug2MassHypo);
495 std::vector<double> massesDx;
496 massesDx.push_back(m_vtx1Daug1MassHypo);
497 massesDx.push_back(m_vtx1Daug2MassHypo);
498 massesDx.push_back(m_vtx1Daug3MassHypo);
499 std::vector<double> massesDm; // Alter the order of masses for D-
500 massesDm.push_back(m_vtx1Daug2MassHypo);
501 massesDm.push_back(m_vtx1Daug1MassHypo);
502 massesDm.push_back(m_vtx1Daug3MassHypo);
503 std::vector<double> Masses;
504 Masses.push_back(m_vtx0Daug1MassHypo);
505 Masses.push_back(m_vtx0Daug2MassHypo);
506 Masses.push_back(m_vtx1MassHypo);
507
508 // Select the J/psi candidates before calling cascade fit
509 std::vector<const xAOD::Vertex*> selectedJpsiCandidates;
510 for(auto vxcItr=jpsiContainer->cbegin(); vxcItr!=jpsiContainer->cend(); ++vxcItr) {
511
512 // Check the passed flag first
513 const xAOD::Vertex* vtx = *vxcItr;
514 SG::AuxElement::Accessor<Char_t> flagAcc1("passed_Jpsi");
515 if(flagAcc1.isAvailable(*vtx)){
516 if(!flagAcc1(*vtx)) continue;
517 }
518
519 // Check J/psi candidate invariant mass and skip if need be
520 double mass_Jpsi = m_V0Tools->invariantMass(*vxcItr, massesJpsi);
521 if (mass_Jpsi < m_jpsiMassLower || mass_Jpsi > m_jpsiMassUpper) {
522 ATH_MSG_DEBUG(" Original Jpsi candidate rejected by the mass cut: mass = "
523 << mass_Jpsi << " != (" << m_jpsiMassLower << ", " << m_jpsiMassUpper << ")" );
524 continue;
525 }
526 selectedJpsiCandidates.push_back(*vxcItr);
527 }
528 if(selectedJpsiCandidates.size()<1) return StatusCode::SUCCESS;
529
530 // Select the D_s+/D+ candidates before calling cascade fit
531 std::vector<const xAOD::Vertex*> selectedDxCandidates;
532 for(auto vxcItr=dxContainer->cbegin(); vxcItr!=dxContainer->cend(); ++vxcItr) {
533
534 // Check the passed flag first
535 const xAOD::Vertex* vtx = *vxcItr;
536 if(abs(m_Dx_pid)==431) { // D_s+/-
537 SG::AuxElement::Accessor<Char_t> flagAcc1("passed_Ds");
538 if(flagAcc1.isAvailable(*vtx)){
539 if(!flagAcc1(*vtx)) continue;
540 }
541 }
542
543 if(abs(m_Dx_pid)==411) { // D+/-
544 SG::AuxElement::Accessor<Char_t> flagAcc1("passed_Dp");
545 SG::AuxElement::Accessor<Char_t> flagAcc2("passed_Dm");
546 bool isDp(true);
547 bool isDm(true);
548 if(flagAcc1.isAvailable(*vtx)){
549 if(!flagAcc1(*vtx)) isDp = false;
550 }
551 if(flagAcc2.isAvailable(*vtx)){
552 if(!flagAcc2(*vtx)) isDm = false;
553 }
554 if(!(isDp||isDm)) continue;
555 }
556
557
558 // Ensure the total charge is correct
559 if(abs((*vxcItr)->trackParticle(0)->charge()+(*vxcItr)->trackParticle(1)->charge()+(*vxcItr)->trackParticle(2)->charge()) != 1){
560 ATH_MSG_DEBUG(" Original D+ candidate rejected by the charge requirement: "
561 << (*vxcItr)->trackParticle(0)->charge() << ", " << (*vxcItr)->trackParticle(1)->charge() << ", " << (*vxcItr)->trackParticle(2)->charge() );
562 continue;
563 }
564
565 // Check D_(s)+/- candidate invariant mass and skip if need be
566 double mass_D;
567 if(abs(m_Dx_pid)==411 && (*vxcItr)->trackParticle(2)->charge()<0) // D-
568 mass_D = m_V0Tools->invariantMass(*vxcItr,massesDm);
569 else // D+, D_s+/-
570 mass_D = m_V0Tools->invariantMass(*vxcItr,massesDx);
571 ATH_MSG_DEBUG("D_(s) mass " << mass_D);
572 if(mass_D < m_DxMassLower || mass_D > m_DxMassUpper) {
573 ATH_MSG_DEBUG(" Original D_(s) candidate rejected by the mass cut: mass = "
574 << mass_D << " != (" << m_DxMassLower << ", " << m_DxMassUpper << ")" );
575 continue;
576 }
577
578 // Add loose cut on K+k- mass for D_s->phi pi
579 if(m_Dx_pid==431){
580 TLorentzVector p4Kp_in, p4Km_in;
581 p4Kp_in.SetPtEtaPhiM( (*vxcItr)->trackParticle(0)->pt(),
582 (*vxcItr)->trackParticle(0)->eta(),
583 (*vxcItr)->trackParticle(0)->phi(), m_vtx1Daug1MassHypo);
584 p4Km_in.SetPtEtaPhiM( (*vxcItr)->trackParticle(1)->pt(),
585 (*vxcItr)->trackParticle(1)->eta(),
586 (*vxcItr)->trackParticle(1)->phi(), m_vtx1Daug2MassHypo);
587 double mass_phi = (p4Kp_in + p4Km_in).M();
588 ATH_MSG_DEBUG("phi mass " << mass_phi);
589 if(mass_phi > 1200) {
590 ATH_MSG_DEBUG(" Original phi candidate rejected by the mass cut: mass = " << mass_phi );
591 continue;
592 }
593 }
594 selectedDxCandidates.push_back(*vxcItr);
595 }
596 if(selectedDxCandidates.size()<1) return StatusCode::SUCCESS;
597
598 // Select J/psi D_(s)+ candidates
599 // Iterate over Jpsi vertices
600 for(auto jpsiItr=selectedJpsiCandidates.cbegin(); jpsiItr!=selectedJpsiCandidates.cend(); ++jpsiItr) {
601
602 size_t jpsiTrkNum = (*jpsiItr)->nTrackParticles();
603 tracksJpsi.clear();
604 for( unsigned int it=0; it<jpsiTrkNum; it++) tracksJpsi.push_back((*jpsiItr)->trackParticle(it));
605
606 if (tracksJpsi.size() != 2 || massesJpsi.size() != 2 ) {
607 ATH_MSG_INFO("problems with Jpsi input");
608 }
609
610 // Iterate over D_(s)+/- vertices
611 for(auto dxItr=selectedDxCandidates.cbegin(); dxItr!=selectedDxCandidates.cend(); ++dxItr) {
612
613 // Check identical tracks in input
614 if(std::find(tracksJpsi.cbegin(), tracksJpsi.cend(), (*dxItr)->trackParticle(0)) != tracksJpsi.cend()) continue;
615 if(std::find(tracksJpsi.cbegin(), tracksJpsi.cend(), (*dxItr)->trackParticle(1)) != tracksJpsi.cend()) continue;
616 if(std::find(tracksJpsi.cbegin(), tracksJpsi.cend(), (*dxItr)->trackParticle(2)) != tracksJpsi.cend()) continue;
617
618 size_t dxTrkNum = (*dxItr)->nTrackParticles();
619 tracksDx.clear();
620 for( unsigned int it=0; it<dxTrkNum; it++) tracksDx.push_back((*dxItr)->trackParticle(it));
621 if (tracksDx.size() != 3 || massesDx.size() != 3 ) {
622 ATH_MSG_INFO("problems with D_(s) input");
623 }
624
625 ATH_MSG_DEBUG("using tracks" << tracksJpsi[0] << ", " << tracksJpsi[1] << ", " << tracksDx[0] << ", " << tracksDx[1] << ", " << tracksDx[2]);
626 tracksBc.clear();
627 for( unsigned int it=0; it<jpsiTrkNum; it++) tracksBc.push_back((*jpsiItr)->trackParticle(it));
628 for( unsigned int it=0; it<dxTrkNum; it++) tracksBc.push_back((*dxItr)->trackParticle(it));
629
630 // Apply the user's settings to the fitter
631 // Reset
632 std::unique_ptr<Trk::IVKalState> state (m_iVertexFitter->makeState());
633 // Robustness
634 int robustness = 0;
635 m_iVertexFitter->setRobustness(robustness, *state);
636 // Build up the topology
637 // Vertex list
638 std::vector<Trk::VertexID> vrtList;
639 // D_(s)+/- vertex
640 Trk::VertexID vID;
641 if (m_constrDx) {
642 if(abs(m_Dx_pid)==411 && (*dxItr)->trackParticle(2)->charge()<0) // D-
643 vID = m_iVertexFitter->startVertex(tracksDx,massesDm,*state,mass_d);
644 else // D+, D_s+/-
645 vID = m_iVertexFitter->startVertex(tracksDx,massesDx,*state,mass_d);
646 } else {
647 if(abs(m_Dx_pid)==411 && (*dxItr)->trackParticle(2)->charge()<0) // D-
648 vID = m_iVertexFitter->startVertex(tracksDx,massesDm,*state);
649 else // D+, D_s+/-
650 vID = m_iVertexFitter->startVertex(tracksDx,massesDx,*state);
651 }
652 vrtList.push_back(vID);
653 // B vertex including Jpsi
654 Trk::VertexID vID2 = m_iVertexFitter->nextVertex(tracksJpsi,massesJpsi,vrtList,*state);
655 if (m_constrJpsi) {
656 std::vector<Trk::VertexID> cnstV;
657 cnstV.clear();
658 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJpsi,cnstV,*state,m_mass_jpsi).isSuccess() ) {
659 ATH_MSG_WARNING("addMassConstraint failed");
660 //return StatusCode::FAILURE;
661 }
662 }
663 // Do the work
664 std::unique_ptr<Trk::VxCascadeInfo> result(m_iVertexFitter->fitCascade(*state));
665
666 if (result != nullptr) {
667 // reset links to original tracks
668 BPhysPVCascadeTools::PrepareVertexLinks(result.get(), trackContainer);
669 ATH_MSG_DEBUG("storing tracks " << ((result->vertices())[0])->trackParticle(0) << ", "
670 << ((result->vertices())[0])->trackParticle(1) << ", "
671 << ((result->vertices())[0])->trackParticle(2) << ", "
672 << ((result->vertices())[1])->trackParticle(0) << ", "
673 << ((result->vertices())[1])->trackParticle(1));
674 // necessary to prevent memory leak
675 result->setSVOwnership(true);
676
677 // Chi2/DOF cut
678 double bChi2DOF = result->fitChi2()/result->nDoF();
679 ATH_MSG_DEBUG("Candidate chi2/DOF is " << bChi2DOF);
680 bool chi2CutPassed = (m_chi2cut <= 0.0 || bChi2DOF < m_chi2cut);
681
682 const std::vector< std::vector<TLorentzVector> > &moms = result->getParticleMoms();
683 double mass = m_CascadeTools->invariantMass(moms[1]);
684 if(chi2CutPassed) {
685 if (mass >= m_MassLower && mass <= m_MassUpper) {
686 cascadeinfoContainer->push_back(result.release());
687 } else {
688 ATH_MSG_DEBUG("Candidate rejected by the mass cut: mass = "
689 << mass << " != (" << m_MassLower << ", " << m_MassUpper << ")" );
690 }
691 }
692 }
693
694 } //Iterate over D_(s)+ vertices
695
696 } //Iterate over Jpsi vertices
697
698 ATH_MSG_DEBUG("cascadeinfoContainer size " << cascadeinfoContainer->size());
699
700 return StatusCode::SUCCESS;
701 }
702
703}
704
705
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_ERROR(x)
#define ATH_MSG_FATAL(x)
#define ATH_MSG_INFO(x)
#define ATH_MSG_WARNING(x)
#define ATH_MSG_DEBUG(x)
#define BPHYS_CHECK(EXP)
Useful CHECK macro.
: B-physics xAOD helpers.
ATLAS-specific HepMC functions.
#define x
const_iterator cbegin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
size_type size() const noexcept
Returns the number of elements in the collection.
const_iterator cend() const noexcept
Return a const_iterator pointing past the end of the collection.
static bool LinkVertices(SG::AuxElement::Decorator< VertexLinkVector > &decor, const std::vector< const xAOD::Vertex * > &vertices, const xAOD::VertexContainer *vertexContainer, const xAOD::Vertex *vert)
static void PrepareVertexLinks(Trk::VxCascadeInfo *result, const xAOD::TrackParticleContainer *importedTrackCollection)
static double getParticleMass(const HepPDT::ParticleDataTable *pdt, int pdg)
static const xAOD::Vertex * FindVertex(const xAOD::VertexContainer *c, const xAOD::Vertex *v)
static void SetVectorInfo(xAOD::BPhysHelper &, const Trk::VxCascadeInfo *)
JpsiPlusDsCascade(const std::string &t, const std::string &n, const IInterface *p)
virtual StatusCode addBranches(const EventContext &ctx) const override
StatusCode performSearch(std::vector< Trk::VxCascadeInfo * > *cascadeinfoContainer, const EventContext &ctx) const
SG::ReadHandleKey< xAOD::EventInfo > m_eventInfo_key
ToolHandle< Analysis::PrimaryVertexRefitter > m_pvRefitter
std::string m_VxPrimaryCandidateName
Name of primary vertex container // FIXME Use Handles.
ToolHandle< Trk::TrkVKalVrtFitter > m_iVertexFitter
std::vector< std::string > m_cascadeOutputsKeys
std::string m_hypoName
name of the mass hypothesis.
ToolHandle< DerivationFramework::CascadeTools > m_CascadeTools
ServiceHandle< IPartPropSvc > m_partPropSvc
virtual StatusCode initialize() override
SG::Decorator< T, ALLOC > Decorator
Definition AuxElement.h:575
SG::Accessor< T, ALLOC > Accessor
Definition AuxElement.h:572
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
bool setMass(const float val)
Set given invariant mass and its error.
bool setMassErr(const float val)
invariant mass error
virtual double phi() const override final
The azimuthal angle ( ) of the particle (has range to .)
virtual double pt() const override final
The transverse momentum ( ) of the particle.
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
float charge() const
Returns the charge.
float z() const
Returns the z position.
const TrackParticle * trackParticle(size_t i) const
Get the pointer to a given track that was used in vertex reco.
float y() const
Returns the y position.
float x() const
Returns the x position.
bool contains(const std::string &s, const std::string &regx)
does a string contain the substring
Definition hcg.cxx:114
THE reconstruction tool.
std::vector< const xAOD::TrackParticle * > TrackBag
ElementLink< xAOD::VertexContainer > VertexLink
std::vector< VertexLink > VertexLinkVector
static const int KPLUS
static const int MUON
static const int DPLUS
static const int BCPLUS
static const int DSPLUS
static const int PIPLUS
static const int JPSI
VertexAuxContainer_v1 VertexAuxContainer
Definition of the current jet auxiliary container.
VertexContainer_v1 VertexContainer
Definition of the current "Vertex container version".
Vertex_v1 Vertex
Define the latest version of the vertex class.
TrackParticleContainer_v1 TrackParticleContainer
Definition of the current "TrackParticle container version".