ATLAS Offline Software
Loading...
Searching...
No Matches
JpsiPlusV0Cascade.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3*/
4
5#include "JpsiPlusV0Cascade.h"
6#include "TrkVertexFitterInterfaces/IVertexFitter.h"
7#include "TrkVKalVrtFitter/TrkVKalVrtFitter.h"
8#include "TrkVertexAnalysisUtils/V0Tools.h"
9#include "DerivationFrameworkBPhys/CascadeTools.h"
10#include "BPhysPVCascadeTools.h"
11#include "xAODTracking/VertexAuxContainer.h"
12#include "xAODEventInfo/EventInfo.h"
13#include "xAODBPhys/BPhysHypoHelper.h"
14#include <algorithm>
15#include "GeneratorModules/GenData.h"
16#include "TrkVKalVrtFitter/VxCascadeInfo.h"
17#include "TruthUtils/HepMCHelpers.h"
18
19namespace DerivationFramework {
20 typedef ElementLink<xAOD::VertexContainer> VertexLink;
21 typedef std::vector<VertexLink> VertexLinkVector;
22 typedef std::vector<const xAOD::TrackParticle*> TrackBag;
23
24
25 JpsiPlusV0Cascade::JpsiPlusV0Cascade(const std::string& t, const std::string& n, const IInterface* p) : base_class(t,n,p)
26 {
27 }
28
29
30 JpsiPlusV0Cascade::~JpsiPlusV0Cascade(){ }
31
32
33 StatusCode JpsiPlusV0Cascade::initialize() {
34 ATH_CHECK(m_eventInfo_key.initialize());
35 ATH_CHECK(m_vertexContainerKey.initialize());
36 ATH_CHECK(m_vertexV0ContainerKey.initialize());
37 ATH_CHECK(m_cascadeOutputsKeys.initialize());
38 ATH_CHECK(m_VxPrimaryCandidateName.initialize());
39 ATH_CHECK(m_jpsiTrackContainerName.initialize());
40 ATH_CHECK(m_v0TrackContainerName.initialize());
41 ATH_CHECK(m_refPVContainerName.initialize());
42
43 // retrieving vertex Fitter
44 ATH_CHECK( m_iVertexFitter.retrieve() );
45 ATH_MSG_DEBUG("Retrieved tool " << m_iVertexFitter);
46
47 // retrieving the V0 tools
48 ATH_CHECK( m_V0Tools.retrieve() );
49 ATH_MSG_INFO("Retrieved tool " << m_V0Tools);
50
51 // retrieving the Cascade tools
52 ATH_CHECK( m_CascadeTools.retrieve() );
53 ATH_MSG_INFO("Retrieved tool " << m_CascadeTools);
54
55 auto gendata = std::make_shared<GenData>();
56 m_mass_electron = gendata->particleMass(MC::ELECTRON).value();
57 m_mass_muon = gendata->particleMass(MC::MUON).value();
58 m_mass_pion = gendata->particleMass(MC::PIPLUS).value();
59 m_mass_proton = gendata->particleMass(MC::PROTON).value();
60 m_mass_lambda = gendata->particleMass(MC::LAMBDA0).value();
61 m_mass_ks = gendata->particleMass(MC::K0S).value();
62 m_mass_jpsi = gendata->particleMass(MC::JPSI).value();
63 m_mass_b0 = gendata->particleMass(MC::B0).value();
64 m_mass_lambdaB = gendata->particleMass(MC::LAMBDAB0).value();
65 ATH_CHECK(m_RelinkContainers.initialize());
66
67 return StatusCode::SUCCESS;
68 }
69
70
71 StatusCode JpsiPlusV0Cascade::addBranches(const EventContext& ctx) const
72 {
73 std::vector<Trk::VxCascadeInfo*> cascadeinfoContainer;
74 constexpr int topoN = 2;
75 std::array<SG::WriteHandle<xAOD::VertexContainer>, topoN> Vtxwritehandles;
76 if (m_cascadeOutputsKeys.size() !=topoN) { ATH_MSG_FATAL("Incorrect number of VtxContainers"); return StatusCode::FAILURE; }
77
78 for (int i =0; i<topoN;i++) {
79 Vtxwritehandles[i] = SG::makeHandle(m_cascadeOutputsKeys[i], ctx);
80 ATH_CHECK( Vtxwritehandles[i].record (std::make_unique<xAOD::VertexContainer>(),
81 std::make_unique<xAOD::VertexAuxContainer>()) );
82 }
83
84 //----------------------------------------------------
85 // retrieve primary vertices
86 //----------------------------------------------------
87 SG::ReadHandle<xAOD::VertexContainer> pvContainer{m_VxPrimaryCandidateName, ctx};
88 if (!pvContainer.isValid()) {
89 ATH_MSG_ERROR("Failed to find xAOD::VertexContainer named " << m_VxPrimaryCandidateName.key() << " in EventStore.");
90 return StatusCode::FAILURE;
91 }
92 ATH_MSG_DEBUG("Found " << m_VxPrimaryCandidateName << " in StoreGate!");
93
94 if (pvContainer->size()==0){
95 ATH_MSG_WARNING("You have no primary vertices: " << pvContainer->size());
96 return StatusCode::RECOVERABLE;
97 }
98 const xAOD::Vertex * primaryVertex = (*pvContainer)[0];
99
100 //----------------------------------------------------
101 // Try to retrieve refitted primary vertices
102 //----------------------------------------------------
103 SG::WriteHandle<xAOD::VertexContainer> refPvContainer;
104 if (m_refitPV) {
105 // refitted PV container does not exist. Create a new one.
106 refPvContainer = SG::makeHandle(m_refPVContainerName, ctx);
107 ATH_CHECK(refPvContainer.record(std::make_unique<xAOD::VertexContainer>(), std::make_unique<xAOD::VertexAuxContainer>()));
108 }
109
110 ATH_CHECK(performSearch(cascadeinfoContainer, ctx));
111
112 SG::ReadHandle<xAOD::EventInfo> evt(m_eventInfo_key, ctx);
113 if(not evt.isValid()) ATH_MSG_ERROR("Cannot Retrieve " << m_eventInfo_key.key() );
114 BPhysPVCascadeTools helper(&(*m_CascadeTools), evt.cptr());
115 helper.SetMinNTracksInPV(m_PV_minNTracks);
116
117 // Decorators for the main vertex: chi2, ndf, pt and pt error, plus the V0 vertex variables
118 SG::AuxElement::Decorator<VertexLinkVector> CascadeLinksDecor("CascadeVertexLinks");
119 SG::AuxElement::Decorator<VertexLinkVector> JpsiLinksDecor("JpsiVertexLinks");
120 SG::AuxElement::Decorator<VertexLinkVector> V0LinksDecor("V0VertexLinks");
121 SG::AuxElement::Decorator<float> chi2_decor("ChiSquared");
122 SG::AuxElement::Decorator<float> ndof_decor("NumberDoF");
123 SG::AuxElement::Decorator<float> Pt_decor("Pt");
124 SG::AuxElement::Decorator<float> PtErr_decor("PtErr");
125 SG::AuxElement::Decorator<float> Mass_svdecor("V0_mass");
126 SG::AuxElement::Decorator<float> MassErr_svdecor("V0_massErr");
127 SG::AuxElement::Decorator<float> Pt_svdecor("V0_Pt");
128 SG::AuxElement::Decorator<float> PtErr_svdecor("V0_PtErr");
129 SG::AuxElement::Decorator<float> Lxy_svdecor("V0_Lxy");
130 SG::AuxElement::Decorator<float> LxyErr_svdecor("V0_LxyErr");
131 SG::AuxElement::Decorator<float> Tau_svdecor("V0_Tau");
132 SG::AuxElement::Decorator<float> TauErr_svdecor("V0_TauErr");
133
134 ATH_MSG_DEBUG("cascadeinfoContainer size " << cascadeinfoContainer.size());
135
136 // Get Jpsi container and identify the input Jpsi
137 SG::ReadHandle<xAOD::VertexContainer>jpsiContainer{m_vertexContainerKey, ctx};
138 if (!jpsiContainer.isValid()) {
139 ATH_MSG_ERROR("Failed to find xAOD::VertexContainer named " << m_vertexContainerKey.key() << " in EventStore.");
140 return StatusCode::FAILURE;
141 }
142 SG::ReadHandle<xAOD::VertexContainer>v0Container{m_vertexV0ContainerKey, ctx};
143 if (!v0Container.isValid()) {
144 ATH_MSG_ERROR("Failed to find xAOD::VertexContainer named " << m_vertexV0ContainerKey.key() << " in EventStore.");
145 return StatusCode::FAILURE;
146 }
147
148 for (Trk::VxCascadeInfo* x : cascadeinfoContainer) {
149 if(x==nullptr) {
150 ATH_MSG_ERROR("cascadeinfoContainer is null");
151 //x is dereferenced if we pass this
152 return StatusCode::FAILURE;
153 }
154
155 // the cascade fitter returns:
156 // std::vector<xAOD::Vertex*>, each xAOD::Vertex contains the refitted track parameters (perigee at the vertex position)
157 // vertices[iv] the links to the original TPs and a covariance of size 3+5*NTRK; the chi2 of the total fit
158 // is split between the cascade vertices as per track contribution
159 // std::vector< std::vector<TLorentzVector> >, each std::vector<TLorentzVector> contains the refitted momenta (TLorentzVector)
160 // momenta[iv][...] of all tracks in the corresponding vertex, including any pseudotracks (from cascade vertices)
161 // originating in this vertex; the masses are as assigned in the cascade fit
162 // std::vector<Amg::MatrixX>, the corresponding covariance matrices in momentum space
163 // covariance[iv]
164 // int nDoF, double Chi2
165 //
166 // the invariant mass, pt, lifetime etc. errors should be calculated using the covariance matrices in momentum space as these
167 // take into account the full track-track and track-vertex correlations
168 //
169 // 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.
170 // The covariance terms between the two vertices are not stored. In momentum space momenta[0] contains the 2 V0 tracks,
171 // their momenta add up to the momentum of the 3rd track in momenta[1], the first two being the Jpsi tracks
172
173 const std::vector<xAOD::Vertex*> &cascadeVertices = x->vertices();
174 if(cascadeVertices.size()!=topoN)
175 ATH_MSG_ERROR("Incorrect number of vertices");
176 if(cascadeVertices[0] == nullptr || cascadeVertices[1] == nullptr) ATH_MSG_ERROR("Error null vertex");
177 // Keep vertices (bear in mind that they come in reverse order!)
178 for(int i =0;i<topoN;i++) Vtxwritehandles[i]->push_back(cascadeVertices[i]);
179
180 x->setSVOwnership(false); // Prevent Container from deleting vertices
181 const auto mainVertex = cascadeVertices[1]; // this is the Bd (Bd, Lambda_b, Lambda_bbar) vertex
182 //const auto v0Vertex = cascadeVertices[0]; // this is the V0 (Kshort, Lambda, Lambdabar) vertex
183 const std::vector< std::vector<TLorentzVector> > &moms = x->getParticleMoms();
184
185 // Set links to cascade vertices
186 std::vector<const xAOD::Vertex*> verticestoLink;
187 verticestoLink.push_back(cascadeVertices[0]);
188 if(!Vtxwritehandles[1].isValid()) ATH_MSG_ERROR("Vtxwritehandles[1] is not valid");
189 if(!BPhysPVCascadeTools::LinkVertices(CascadeLinksDecor, verticestoLink, Vtxwritehandles[0].cptr(), cascadeVertices[1]))
190 ATH_MSG_ERROR("Error decorating with cascade vertices");
191
192 // Identify the input Jpsi
193 const xAOD::Vertex* jpsiVertex = BPhysPVCascadeTools::FindVertex<2>(jpsiContainer.cptr(), cascadeVertices[1]);
194 ATH_MSG_DEBUG("1 pt Jpsi tracks " << cascadeVertices[1]->trackParticle(0)->pt() << ", " << cascadeVertices[1]->trackParticle(1)->pt());
195 if (jpsiVertex) ATH_MSG_DEBUG("2 pt Jpsi tracks " << jpsiVertex->trackParticle(0)->pt() << ", " << jpsiVertex->trackParticle(1)->pt());
196
197 // Identify the input V0
198 const xAOD::Vertex* v0Vertex = BPhysPVCascadeTools::FindVertex<2>(v0Container.cptr(), cascadeVertices[0]);;
199 ATH_MSG_DEBUG("1 pt V0 tracks " << cascadeVertices[0]->trackParticle(0)->pt() << ", " << cascadeVertices[0]->trackParticle(1)->pt());
200 if (v0Vertex) ATH_MSG_DEBUG("2 pt V0 tracks " << v0Vertex->trackParticle(0)->pt() << ", " << v0Vertex->trackParticle(1)->pt());
201
202 // Set links to input vertices
203 std::vector<const xAOD::Vertex*> jpsiVerticestoLink;
204 if (jpsiVertex) jpsiVerticestoLink.push_back(jpsiVertex);
205 else ATH_MSG_WARNING("Could not find linking Jpsi");
206 if(!BPhysPVCascadeTools::LinkVertices(JpsiLinksDecor, jpsiVerticestoLink, jpsiContainer.cptr(), cascadeVertices[1]))
207 ATH_MSG_ERROR("Error decorating with Jpsi vertices");
208
209 std::vector<const xAOD::Vertex*> v0VerticestoLink;
210 if (v0Vertex) v0VerticestoLink.push_back(v0Vertex);
211 else ATH_MSG_WARNING("Could not find linking V0");
212 if(!BPhysPVCascadeTools::LinkVertices(V0LinksDecor, v0VerticestoLink, v0Container.cptr(), cascadeVertices[1]))
213 ATH_MSG_ERROR("Error decorating with V0 vertices");
214
215 double mass_v0 = m_mass_ks;
216 double mass_b = m_mass_b0;
217 double mass_track = MC::isElectron(m_jpsi_trk_pdg.value()) ? m_mass_electron : m_mass_muon;
218 std::vector<double> massesJpsi(2, mass_track);
219 std::vector<double> massesV0;
220 std::vector<double> Masses(2, mass_track);
221 if (m_v0_pid == 310) {
222 massesV0.push_back(m_mass_pion);
223 massesV0.push_back(m_mass_pion);
224 Masses.push_back(m_mass_ks);
225 } else if (m_v0_pid == 3122) {
226 massesV0.push_back(m_mass_proton);
227 massesV0.push_back(m_mass_pion);
228 Masses.push_back(m_mass_lambda);
229 mass_v0 = m_mass_lambda;
230 mass_b = m_mass_lambdaB;
231 } else if (m_v0_pid == -3122) {
232 massesV0.push_back(m_mass_pion);
233 massesV0.push_back(m_mass_proton);
234 Masses.push_back(m_mass_lambda);
235 mass_v0 = m_mass_lambda;
236 mass_b = m_mass_lambdaB;
237 }
238
239 // loop over candidates -- Don't apply PV_minNTracks requirement here
240 // because it may result in exclusion of the high-pt PV.
241 // get good PVs
242
243 xAOD::BPhysHypoHelper vtx(m_hypoName, mainVertex);
244
245 BPhysPVCascadeTools::SetVectorInfo(vtx, x);
246
247
248 // Decorate main vertex
249 //
250 // 1.a) mass, mass error
251 BPHYS_CHECK( vtx.setMass(m_CascadeTools->invariantMass(moms[1])) );
252 BPHYS_CHECK( vtx.setMassErr(m_CascadeTools->invariantMassError(moms[1],x->getCovariance()[1])) );
253 // 1.b) pt and pT error (the default pt of mainVertex is != the pt of the full cascade fit!)
254 Pt_decor(*mainVertex) = m_CascadeTools->pT(moms[1]);
255 PtErr_decor(*mainVertex) = m_CascadeTools->pTError(moms[1],x->getCovariance()[1]);
256 // 1.c) chi2 and ndof (the default chi2 of mainVertex is != the chi2 of the full cascade fit!)
257 chi2_decor(*mainVertex) = x->fitChi2();
258 ndof_decor(*mainVertex) = x->nDoF();
259
260 xAOD::VertexContainer *refPVContainer{};
261 if (m_refitPV) { refPVContainer = refPvContainer.ptr(); }
262 ATH_CHECK(helper.FillCandwithRefittedVertices(m_refitPV, pvContainer.cptr(),
263 refPVContainer, &(*m_pvRefitter), m_PV_max, m_DoVertexType, x, 1, mass_b, vtx));
264
265
266 // 4) decorate the main vertex with V0 vertex mass, pt, lifetime and lxy values (plus errors)
267 // V0 points to the main vertex, so lifetime and lxy are w.r.t the main vertex
268 Mass_svdecor(*mainVertex) = m_CascadeTools->invariantMass(moms[0]);
269 MassErr_svdecor(*mainVertex) = m_CascadeTools->invariantMassError(moms[0],x->getCovariance()[0]);
270 Pt_svdecor(*mainVertex) = m_CascadeTools->pT(moms[0]);
271 PtErr_svdecor(*mainVertex) = m_CascadeTools->pTError(moms[0],x->getCovariance()[0]);
272 Lxy_svdecor(*mainVertex) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
273 LxyErr_svdecor(*mainVertex) = m_CascadeTools->lxyError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
274 Tau_svdecor(*mainVertex) = m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1]);
275 TauErr_svdecor(*mainVertex) = m_CascadeTools->tauError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
276
277 // Some checks in DEBUG mode
278 ATH_MSG_DEBUG("chi2 " << x->fitChi2()
279 << " chi2_1 " << m_V0Tools->chisq(cascadeVertices[0])
280 << " chi2_2 " << m_V0Tools->chisq(cascadeVertices[1])
281 << " vprob " << m_CascadeTools->vertexProbability(x->nDoF(),x->fitChi2()));
282 ATH_MSG_DEBUG("ndf " << x->nDoF() << " ndf_1 " << m_V0Tools->ndof(cascadeVertices[0]) << " ndf_2 " << m_V0Tools->ndof(cascadeVertices[1]));
283 ATH_MSG_DEBUG("V0Tools mass_v0 " << m_V0Tools->invariantMass(cascadeVertices[0],massesV0)
284 << " error " << m_V0Tools->invariantMassError(cascadeVertices[0],massesV0)
285 << " mass_J " << m_V0Tools->invariantMass(cascadeVertices[1],massesJpsi)
286 << " error " << m_V0Tools->invariantMassError(cascadeVertices[1],massesJpsi));
287 // masses and errors, using track masses assigned in the fit
288 double Mass_B = m_CascadeTools->invariantMass(moms[1]);
289 double Mass_V0 = m_CascadeTools->invariantMass(moms[0]);
290 double Mass_B_err = m_CascadeTools->invariantMassError(moms[1],x->getCovariance()[1]);
291 double Mass_V0_err = m_CascadeTools->invariantMassError(moms[0],x->getCovariance()[0]);
292 ATH_MSG_DEBUG("Mass_B " << Mass_B << " Mass_V0 " << Mass_V0);
293 ATH_MSG_DEBUG("Mass_B_err " << Mass_B_err << " Mass_V0_err " << Mass_V0_err);
294 double mprob_B = m_CascadeTools->massProbability(mass_b,Mass_B,Mass_B_err);
295 double mprob_V0 = m_CascadeTools->massProbability(mass_v0,Mass_V0,Mass_V0_err);
296 ATH_MSG_DEBUG("mprob_B " << mprob_B << " mprob_V0 " << mprob_V0);
297 // masses and errors, assigning user defined track masses
298 ATH_MSG_DEBUG("Mass_b " << m_CascadeTools->invariantMass(moms[1],Masses)
299 << " Mass_v0 " << m_CascadeTools->invariantMass(moms[0],massesV0));
300 ATH_MSG_DEBUG("Mass_b_err " << m_CascadeTools->invariantMassError(moms[1],x->getCovariance()[1],Masses)
301 << " Mass_v0_err " << m_CascadeTools->invariantMassError(moms[0],x->getCovariance()[0],massesV0));
302 ATH_MSG_DEBUG("pt_b " << m_CascadeTools->pT(moms[1])
303 << " pt_v " << m_CascadeTools->pT(moms[0])
304 << " pt_v0 " << m_V0Tools->pT(cascadeVertices[0]));
305 ATH_MSG_DEBUG("ptErr_b " << m_CascadeTools->pTError(moms[1],x->getCovariance()[1])
306 << " ptErr_v " << m_CascadeTools->pTError(moms[0],x->getCovariance()[0])
307 << " ptErr_v0 " << m_V0Tools->pTError(cascadeVertices[0]));
308 ATH_MSG_DEBUG("lxy_B " << m_V0Tools->lxy(cascadeVertices[1],primaryVertex) << " lxy_V " << m_V0Tools->lxy(cascadeVertices[0],cascadeVertices[1]));
309 ATH_MSG_DEBUG("lxy_b " << m_CascadeTools->lxy(moms[1],cascadeVertices[1],primaryVertex) << " lxy_v " << m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]));
310 ATH_MSG_DEBUG("lxyErr_b " << m_CascadeTools->lxyError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
311 << " lxyErr_v " << m_CascadeTools->lxyError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
312 << " lxyErr_v0 " << m_V0Tools->lxyError(cascadeVertices[0],cascadeVertices[1]));
313 ATH_MSG_DEBUG("tau_B " << m_CascadeTools->tau(moms[1],cascadeVertices[1],primaryVertex,mass_b)
314 << " tau_v0 " << m_V0Tools->tau(cascadeVertices[0],cascadeVertices[1],massesV0));
315 ATH_MSG_DEBUG("tau_b " << m_CascadeTools->tau(moms[1],cascadeVertices[1],primaryVertex)
316 << " tau_v " << m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1])
317 << " tau_V " << m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1],mass_v0));
318 ATH_MSG_DEBUG("tauErr_b " << m_CascadeTools->tauError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
319 << " tauErr_v " << m_CascadeTools->tauError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
320 << " tauErr_v0 " << m_V0Tools->tauError(cascadeVertices[0],cascadeVertices[1],massesV0));
321 ATH_MSG_DEBUG("TauErr_b " << m_CascadeTools->tauError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex,mass_b)
322 << " TauErr_v " << m_CascadeTools->tauError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1],mass_v0)
323 << " TauErr_v0 " << m_V0Tools->tauError(cascadeVertices[0],cascadeVertices[1],massesV0,mass_v0));
324
325 ATH_MSG_DEBUG("CascadeTools main vert wrt PV " << " CascadeTools SV " << " V0Tools SV");
326 ATH_MSG_DEBUG("a0z " << m_CascadeTools->a0z(moms[1],cascadeVertices[1],primaryVertex)
327 << ", " << m_CascadeTools->a0z(moms[0],cascadeVertices[0],cascadeVertices[1])
328 << ", " << m_V0Tools->a0z(cascadeVertices[0],cascadeVertices[1]));
329 ATH_MSG_DEBUG("a0zErr " << m_CascadeTools->a0zError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
330 << ", " << m_CascadeTools->a0zError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
331 << ", " << m_V0Tools->a0zError(cascadeVertices[0],cascadeVertices[1]));
332 ATH_MSG_DEBUG("a0xy " << m_CascadeTools->a0xy(moms[1],cascadeVertices[1],primaryVertex)
333 << ", " << m_CascadeTools->a0xy(moms[0],cascadeVertices[0],cascadeVertices[1])
334 << ", " << m_V0Tools->a0xy(cascadeVertices[0],cascadeVertices[1]));
335 ATH_MSG_DEBUG("a0xyErr " << m_CascadeTools->a0xyError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
336 << ", " << m_CascadeTools->a0xyError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
337 << ", " << m_V0Tools->a0xyError(cascadeVertices[0],cascadeVertices[1]));
338 ATH_MSG_DEBUG("a0 " << m_CascadeTools->a0(moms[1],cascadeVertices[1],primaryVertex)
339 << ", " << m_CascadeTools->a0(moms[0],cascadeVertices[0],cascadeVertices[1])
340 << ", " << m_V0Tools->a0(cascadeVertices[0],cascadeVertices[1]));
341 ATH_MSG_DEBUG("a0Err " << m_CascadeTools->a0Error(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
342 << ", " << m_CascadeTools->a0Error(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
343 << ", " << m_V0Tools->a0Error(cascadeVertices[0],cascadeVertices[1]));
344 ATH_MSG_DEBUG("x0 " << m_V0Tools->vtx(cascadeVertices[0]).x() << " y0 " << m_V0Tools->vtx(cascadeVertices[0]).y() << " z0 " << m_V0Tools->vtx(cascadeVertices[0]).z());
345 ATH_MSG_DEBUG("x1 " << m_V0Tools->vtx(cascadeVertices[1]).x() << " y1 " << m_V0Tools->vtx(cascadeVertices[1]).y() << " z1 " << m_V0Tools->vtx(cascadeVertices[1]).z());
346 ATH_MSG_DEBUG("X0 " << primaryVertex->x() << " Y0 " << primaryVertex->y() << " Z0 " << primaryVertex->z());
347 ATH_MSG_DEBUG("rxy0 " << m_V0Tools->rxy(cascadeVertices[0]) << " rxyErr0 " << m_V0Tools->rxyError(cascadeVertices[0]));
348 ATH_MSG_DEBUG("rxy1 " << m_V0Tools->rxy(cascadeVertices[1]) << " rxyErr1 " << m_V0Tools->rxyError(cascadeVertices[1]));
349 ATH_MSG_DEBUG("Rxy0 wrt PV " << m_V0Tools->rxy(cascadeVertices[0],primaryVertex) << " RxyErr0 wrt PV " << m_V0Tools->rxyError(cascadeVertices[0],primaryVertex));
350 ATH_MSG_DEBUG("Rxy1 wrt PV " << m_V0Tools->rxy(cascadeVertices[1],primaryVertex) << " RxyErr1 wrt PV " << m_V0Tools->rxyError(cascadeVertices[1],primaryVertex));
351 ATH_MSG_DEBUG("number of covariance matrices " << (x->getCovariance()).size());
352 } // loop over cascadeinfoContainer
353
354 // Deleting cascadeinfo since this won't be stored.
355 // Vertices have been kept in m_cascadeOutputs and should be owned by their container
356 for (auto x : cascadeinfoContainer) delete x;
357
358 return StatusCode::SUCCESS;
359 }
360
361
362 StatusCode JpsiPlusV0Cascade::performSearch(std::vector<Trk::VxCascadeInfo*>& cascadeinfoContainer, const EventContext& ctx) const
363 {
364 ATH_MSG_DEBUG( "JpsiPlusV0Cascade::performSearch" );
365
366 // Get TrackParticle containers (for setting links to the original tracks)
367 SG::ReadHandle<xAOD::TrackParticleContainer> jpsiTrackContainer{m_jpsiTrackContainerName, ctx};
368 if (!jpsiTrackContainer.isValid()) {
369 ATH_MSG_ERROR("Failed to find xAOD::TrackParticleContainer named " << m_jpsiTrackContainerName.key() << " in EventStore.");
370 return StatusCode::FAILURE;
371 }
372 SG::ReadHandle<xAOD::TrackParticleContainer> v0TrackContainer{m_v0TrackContainerName, ctx};
373 if (!v0TrackContainer.isValid()) {
374 ATH_MSG_ERROR("Failed to find xAOD::TrackParticleContainer named " << m_v0TrackContainerName.key() << " in EventStore.");
375 return StatusCode::FAILURE;
376 }
377
378 // Get Jpsi container - TODO might be cleaner to retrieve these once and pass to this method?
379 SG::ReadHandle<xAOD::VertexContainer>jpsiContainer{m_vertexContainerKey, ctx};
380 if (!jpsiContainer.isValid()) {
381 ATH_MSG_ERROR("Failed to find xAOD::VertexContainer named " << m_vertexContainerKey.key() << " in EventStore.");
382 return StatusCode::FAILURE;
383 }
384 SG::ReadHandle<xAOD::VertexContainer>v0Container{m_vertexV0ContainerKey, ctx};
385 if (!v0Container.isValid()) {
386 ATH_MSG_ERROR("Failed to find xAOD::VertexContainer named " << m_vertexV0ContainerKey.key() << " in EventStore.");
387 return StatusCode::FAILURE;
388 }
389
390 double mass_v0 = m_mass_ks;
391 double mass_tracks = MC::isElectron(m_jpsi_trk_pdg.value()) ? m_mass_electron : m_mass_muon;
392 std::vector<const xAOD::TrackParticle*> tracksJpsi;
393 std::vector<const xAOD::TrackParticle*> tracksV0;
394 std::vector<double> massesJpsi(2, mass_tracks);
395 std::vector<double> massesV0;
396 std::vector<double> Masses(2, mass_tracks);
397 if (m_v0_pid == 310) {
398 massesV0.push_back(m_mass_pion);
399 massesV0.push_back(m_mass_pion);
400 Masses.push_back(m_mass_ks);
401 } else if (m_v0_pid == 3122) {
402 massesV0.push_back(m_mass_proton);
403 massesV0.push_back(m_mass_pion);
404 mass_v0 = m_mass_lambda;
405 Masses.push_back(m_mass_lambda);
406 } else if (m_v0_pid == -3122) {
407 massesV0.push_back(m_mass_pion);
408 massesV0.push_back(m_mass_proton);
409 mass_v0 = m_mass_lambda;
410 Masses.push_back(m_mass_lambda);
411 }
412 std::vector<const xAOD::TrackParticleContainer*> trackCols;
413 for(const auto &str : m_RelinkContainers){
414 SG::ReadHandle<xAOD::TrackParticleContainer> handle(str,ctx);
415 trackCols.push_back(handle.cptr());
416 }
417
418
419 for(auto jpsi : *jpsiContainer) { //Iterate over Jpsi vertices
420
421 size_t jpsiTrkNum = jpsi->nTrackParticles();
422 tracksJpsi.clear();
423 for( unsigned int it=0; it<jpsiTrkNum; it++) tracksJpsi.push_back(jpsi->trackParticle(it));
424
425 if (tracksJpsi.size() != 2 || massesJpsi.size() != 2 ) {
426 ATH_MSG_INFO("problems with Jpsi input");
427 }
428 double mass_Jpsi = m_V0Tools->invariantMass(jpsi,massesJpsi);
429 ATH_MSG_DEBUG("Jpsi mass " << mass_Jpsi);
430 if (mass_Jpsi < m_jpsiMassLower || mass_Jpsi > m_jpsiMassUpper) {
431 ATH_MSG_DEBUG(" Original Jpsi candidate rejected by the mass cut: mass = "
432 << mass_Jpsi << " != (" << m_jpsiMassLower << ", " << m_jpsiMassUpper << ")" );
433 continue;
434 }
435
436 for(auto v0 : *v0Container) { //Iterate over V0 vertices
437
438 size_t v0TrkNum = v0->nTrackParticles();
439 tracksV0.clear();
440 for( unsigned int it=0; it<v0TrkNum; it++) tracksV0.push_back(v0->trackParticle(it));
441 if (tracksV0.size() != 2 || massesV0.size() != 2 ) {
442 ATH_MSG_INFO("problems with V0 input");
443 }
444 double mass_V0 = m_V0Tools->invariantMass(v0,massesV0);
445 ATH_MSG_DEBUG("V0 mass " << mass_V0);
446 if (mass_V0 < m_V0MassLower || mass_V0 > m_V0MassUpper) {
447 ATH_MSG_DEBUG(" Original V0 candidate rejected by the mass cut: mass = "
448 << mass_V0 << " != (" << m_V0MassLower << ", " << m_V0MassUpper << ")" );
449 continue;
450 }
451 ATH_MSG_DEBUG("using tracks" << tracksJpsi[0] << ", " << tracksJpsi[1] << ", " << tracksV0[0] << ", " << tracksV0[1]);
452 if(!BPhysPVCascadeTools::uniqueCollection(tracksJpsi, tracksV0)) continue;
453
454 // Apply the user's settings to the fitter
455 // Reset
456 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
457 // Robustness
458 int robustness = 0;
459 m_iVertexFitter->setRobustness(robustness, *state);
460 // Build up the topology
461 // Vertex list
462 std::vector<Trk::VertexID> vrtList;
463 // V0 vertex
464 Trk::VertexID vID;
465 if (m_constrV0) {
466 vID = m_iVertexFitter->startVertex(tracksV0,massesV0,*state, mass_v0);
467 } else {
468 vID = m_iVertexFitter->startVertex(tracksV0,massesV0, *state);
469 }
470 vrtList.push_back(vID);
471 // B vertex including Jpsi
472 Trk::VertexID vID2 = m_iVertexFitter->nextVertex(tracksJpsi,massesJpsi,vrtList, *state);
473 if (m_constrJpsi) {
474 std::vector<Trk::VertexID> cnstV;
475 cnstV.clear();
476 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJpsi,cnstV,*state, m_mass_jpsi).isSuccess() ) {
477 ATH_MSG_WARNING("addMassConstraint failed");
478 //return StatusCode::FAILURE;
479 }
480 }
481 // Do the work
482 std::unique_ptr<Trk::VxCascadeInfo> result(m_iVertexFitter->fitCascade(*state));
483
484 if (result) {
485 // reset links to original tracks
486 if(trackCols.empty()) BPhysPVCascadeTools::PrepareVertexLinks(result.get(), v0TrackContainer.cptr());
487 else BPhysPVCascadeTools::PrepareVertexLinks(result.get(), trackCols);
488
489 ATH_MSG_DEBUG("storing tracks " << ((result->vertices())[0])->trackParticle(0) << ", "
490 << ((result->vertices())[0])->trackParticle(1) << ", "
491 << ((result->vertices())[1])->trackParticle(0) << ", "
492 << ((result->vertices())[1])->trackParticle(1));
493
494 // necessary to prevent memory leak
495 result->setSVOwnership(true);
496 const std::vector< std::vector<TLorentzVector> > &moms = result->getParticleMoms();
497 if(moms.size() < 2){
498 ATH_MSG_FATAL("Incorrect size " << __FILE__ << __LINE__ );
499 return StatusCode::FAILURE;
500 }
501 double mass = m_CascadeTools->invariantMass(moms[1]);
502 if (mass >= m_MassLower && mass <= m_MassUpper) {
503
504 cascadeinfoContainer.push_back(result.release());
505 } else {
506 ATH_MSG_DEBUG("Candidate rejected by the mass cut: mass = "
507 << mass << " != (" << m_MassLower << ", " << m_MassUpper << ")" );
508 }
509 }
510
511 } //Iterate over V0 vertices
512
513 } //Iterate over Jpsi vertices
514
515 ATH_MSG_DEBUG("cascadeinfoContainer size " << cascadeinfoContainer.size());
516
517 return StatusCode::SUCCESS;
518 }
519
520
521
522}
ATH_CHECK
#define ATH_CHECK
Evaluate an expression and check for errors.
Definition AthCheckMacros.h:40
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition AthMsgStreamMacros.h:33
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition AthMsgStreamMacros.h:34
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition AthMsgStreamMacros.h:31
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition AthMsgStreamMacros.h:32
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
BPHYS_CHECK
#define BPHYS_CHECK(EXP)
Useful CHECK macro.
Definition BPhysHelper.h:738
BPhysHypoHelper.h
: B-physics xAOD helpers.
HepMCHelpers.h
ATLAS-specific HepMC functions.
x
#define x
DerivationFramework::BPhysPVCascadeTools::uniqueCollection
static bool uniqueCollection(const std::vector< const xAOD::TrackParticle * > &)
Definition BPhysPVCascadeTools.cxx:442
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:459
DerivationFramework::BPhysPVCascadeTools::PrepareVertexLinks
static void PrepareVertexLinks(Trk::VxCascadeInfo *result, const xAOD::TrackParticleContainer *importedTrackCollection)
Definition BPhysPVCascadeTools.cxx:203
DerivationFramework::BPhysPVCascadeTools::FindVertex
static const xAOD::Vertex * FindVertex(const xAOD::VertexContainer *c, const xAOD::Vertex *v)
Definition BPhysPVCascadeTools.h:140
DerivationFramework::BPhysPVCascadeTools::SetVectorInfo
static void SetVectorInfo(xAOD::BPhysHelper &, const Trk::VxCascadeInfo *)
Definition BPhysPVCascadeTools.cxx:423
DerivationFramework::JpsiPlusV0Cascade::m_iVertexFitter
PublicToolHandle< Trk::TrkVKalVrtFitter > m_iVertexFitter
Definition JpsiPlusV0Cascade.h:68
DerivationFramework::JpsiPlusV0Cascade::m_pvRefitter
PublicToolHandle< Analysis::PrimaryVertexRefitter > m_pvRefitter
Definition JpsiPlusV0Cascade.h:69
DerivationFramework::JpsiPlusV0Cascade::m_V0MassUpper
Gaudi::Property< double > m_V0MassUpper
Definition JpsiPlusV0Cascade.h:61
DerivationFramework::JpsiPlusV0Cascade::m_vertexContainerKey
SG::ReadHandleKey< xAOD::VertexContainer > m_vertexContainerKey
Definition JpsiPlusV0Cascade.h:50
DerivationFramework::JpsiPlusV0Cascade::m_VxPrimaryCandidateName
SG::ReadHandleKey< xAOD::VertexContainer > m_VxPrimaryCandidateName
Name of primary vertex container.
Definition JpsiPlusV0Cascade.h:53
DerivationFramework::JpsiPlusV0Cascade::m_CascadeTools
PublicToolHandle< DerivationFramework::CascadeTools > m_CascadeTools
Definition JpsiPlusV0Cascade.h:71
DerivationFramework::JpsiPlusV0Cascade::JpsiPlusV0Cascade
JpsiPlusV0Cascade(const std::string &t, const std::string &n, const IInterface *p)
Definition JpsiPlusV0Cascade.cxx:25
DerivationFramework::JpsiPlusV0Cascade::m_jpsiTrackContainerName
SG::ReadHandleKey< xAOD::TrackParticleContainer > m_jpsiTrackContainerName
Definition JpsiPlusV0Cascade.h:76
DerivationFramework::JpsiPlusV0Cascade::addBranches
virtual StatusCode addBranches(const EventContext &ctx) const override
Definition JpsiPlusV0Cascade.cxx:71
DerivationFramework::JpsiPlusV0Cascade::m_V0MassLower
Gaudi::Property< double > m_V0MassLower
Definition JpsiPlusV0Cascade.h:60
DerivationFramework::JpsiPlusV0Cascade::m_RelinkContainers
SG::ReadHandleKeyArray< xAOD::TrackParticleContainer > m_RelinkContainers
Definition JpsiPlusV0Cascade.h:54
DerivationFramework::JpsiPlusV0Cascade::m_V0Tools
PublicToolHandle< Trk::V0Tools > m_V0Tools
Definition JpsiPlusV0Cascade.h:70
DerivationFramework::JpsiPlusV0Cascade::performSearch
StatusCode performSearch(std::vector< Trk::VxCascadeInfo * > &cascadeinfoContainer, const EventContext &ctx) const
Definition JpsiPlusV0Cascade.cxx:362
DerivationFramework::JpsiPlusV0Cascade::m_vertexV0ContainerKey
SG::ReadHandleKey< xAOD::VertexContainer > m_vertexV0ContainerKey
Definition JpsiPlusV0Cascade.h:51
DerivationFramework::JpsiPlusV0Cascade::m_eventInfo_key
SG::ReadHandleKey< xAOD::EventInfo > m_eventInfo_key
Definition JpsiPlusV0Cascade.h:49
DerivationFramework::JpsiPlusV0Cascade::m_refPVContainerName
SG::WriteHandleKey< xAOD::VertexContainer > m_refPVContainerName
Definition JpsiPlusV0Cascade.h:75
DerivationFramework::JpsiPlusV0Cascade::initialize
virtual StatusCode initialize() override
Definition JpsiPlusV0Cascade.cxx:33
DerivationFramework::JpsiPlusV0Cascade::m_PV_minNTracks
Gaudi::Property< size_t > m_PV_minNTracks
Definition JpsiPlusV0Cascade.h:83
DerivationFramework::JpsiPlusV0Cascade::m_jpsiMassUpper
Gaudi::Property< double > m_jpsiMassUpper
Definition JpsiPlusV0Cascade.h:59
DerivationFramework::JpsiPlusV0Cascade::m_hypoName
Gaudi::Property< std::string > m_hypoName
name of the mass hypothesis.
Definition JpsiPlusV0Cascade.h:78
DerivationFramework::JpsiPlusV0Cascade::m_jpsiMassLower
Gaudi::Property< double > m_jpsiMassLower
Definition JpsiPlusV0Cascade.h:58
DerivationFramework::JpsiPlusV0Cascade::m_v0TrackContainerName
SG::ReadHandleKey< xAOD::TrackParticleContainer > m_v0TrackContainerName
Definition JpsiPlusV0Cascade.h:77
DerivationFramework::JpsiPlusV0Cascade::m_cascadeOutputsKeys
SG::WriteHandleKeyArray< xAOD::VertexContainer > m_cascadeOutputsKeys
Definition JpsiPlusV0Cascade.h:56
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
SG::ReadHandle::cptr
const_pointer_type cptr()
Dereference the pointer.
SG::WriteHandle::record
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
SG::WriteHandle::ptr
pointer_type ptr()
Dereference the pointer.
xAOD::BPhysHypoHelper::setMass
bool setMass(const float val)
Set given invariant mass and its error.
Definition BPhysHypoHelper.cxx:52
xAOD::BPhysHypoHelper::setMassErr
bool setMassErr(const float val)
invariant mass error
Definition BPhysHypoHelper.cxx:57
xAOD::TrackParticle_v1::pt
virtual double pt() const override final
The transverse momentum ( ) of the particle.
Definition TrackParticle_v1.cxx:75
xAOD::Vertex_v1::z
float z() const
Returns the z position.
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:271
xAOD::Vertex_v1::y
float y() const
Returns the y position.
xAOD::Vertex_v1::x
float x() const
Returns the x position.
DerivationFramework
THE reconstruction tool.
Definition CascadeTools.h:16
DerivationFramework::TrackBag
std::vector< const xAOD::TrackParticle * > TrackBag
Definition JpsiPlusDpstCascade.cxx:25
DerivationFramework::VertexLink
ElementLink< xAOD::VertexContainer > VertexLink
Definition JpsiPlusDpstCascade.cxx:23
DerivationFramework::VertexLinkVector
std::vector< VertexLink > VertexLinkVector
Definition JpsiPlusDpstCascade.cxx:24
MC::MUON
static const int MUON
Definition HepMCHelpers.h:78
MC::isElectron
bool isElectron(const T &p)
Definition HepMCHelpers.h:203
MC::ELECTRON
static const int ELECTRON
Definition HepMCHelpers.h:75
MC::K0S
static const int K0S
Definition HepMCHelpers.h:115
MC::LAMBDAB0
static const int LAMBDAB0
Definition HepMCHelpers.h:129
MC::PIPLUS
static const int PIPLUS
Definition HepMCHelpers.h:110
MC::JPSI
static const int JPSI
Definition HepMCHelpers.h:122
MC::B0
static const int B0
Definition HepMCHelpers.h:123
MC::LAMBDA0
static const int LAMBDA0
Definition HepMCHelpers.h:127
MC::PROTON
static const int PROTON
Definition HepMCHelpers.h:125
SG::makeHandle
SG::ReadCondHandle< T > makeHandle(const SG::ReadCondHandleKey< T > &key, const EventContext &ctx=Gaudi::Hive::currentContext())
Definition ReadCondHandle.h:269
xAOD::VertexContainer
VertexContainer_v1 VertexContainer
Definition of the current "Vertex container version".
Definition VertexContainer.h:14
xAOD::Vertex
Vertex_v1 Vertex
Define the latest version of the vertex class.
Definition Event/xAOD/xAODTracking/xAODTracking/Vertex.h:16
Generated on for ATLAS Offline Software by doxygen 1.16.1