ATLAS Offline Software
Loading...
Searching...
No Matches
PsiPlusPsiCascade.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3 Contact: Xin Chen <xin.chen@cern.ch>
4*/
5#include "PsiPlusPsiCascade.h"
11#include "BPhysPVCascadeTools.h"
16#include <algorithm>
17#include <functional>
18
19namespace DerivationFramework {
20 typedef ElementLink<xAOD::VertexContainer> VertexLink;
21 typedef std::vector<VertexLink> VertexLinkVector;
22
24 // retrieving vertex Fitter
25 ATH_CHECK( m_iVertexFitter.retrieve() );
26
27 // retrieve PV refitter
28 ATH_CHECK( m_pvRefitter.retrieve() );
29
30 // retrieving the V0 tools
31 ATH_CHECK( m_V0Tools.retrieve() );
32
33 // retrieving the Cascade tools
34 ATH_CHECK( m_CascadeTools.retrieve() );
35
36 ATH_CHECK( m_vertexPsi1ContainerKey.initialize() );
37 ATH_CHECK( m_vertexPsi2ContainerKey.initialize() );
38 ATH_CHECK( m_VxPrimaryCandidateName.initialize() );
39 ATH_CHECK( m_trackContainerName.initialize() );
40 ATH_CHECK( m_refPVContainerName.initialize() );
41 ATH_CHECK( m_cascadeOutputsKeys.initialize() );
42 ATH_CHECK( m_eventInfo_key.initialize() );
43
44 auto gendata = std::make_shared<GenData>();
45 if(m_massJpsi1 < 0.) m_massJpsi1 = gendata->particleMass(MC::JPSI).value();
46 if(m_massJpsi2 < 0.) m_massJpsi2 = gendata->particleMass(MC::JPSI).value();
47 if(m_massPsi1 < 0.) m_massPsi1 = gendata->particleMass(MC::PSI2S).value();
48 if(m_massPsi2 < 0.) m_massPsi2 = gendata->particleMass(MC::PSI2S).value();
49 if(m_vtx1Daug1MassHypo < 0.) m_vtx1Daug1MassHypo = gendata->particleMass(MC::MUON).value();
50 if(m_vtx1Daug2MassHypo < 0.) m_vtx1Daug2MassHypo = gendata->particleMass(MC::MUON).value();
51 if(m_vtx1Daug3MassHypo < 0.) m_vtx1Daug3MassHypo = gendata->particleMass(MC::PIPLUS).value();
52 if(m_vtx1Daug4MassHypo < 0.) m_vtx1Daug4MassHypo = gendata->particleMass(MC::PIPLUS).value();
53 if(m_vtx2Daug1MassHypo < 0.) m_vtx2Daug1MassHypo = gendata->particleMass(MC::MUON).value();
54 if(m_vtx2Daug2MassHypo < 0.) m_vtx2Daug2MassHypo = gendata->particleMass(MC::MUON).value();
55 if(m_vtx2Daug3MassHypo < 0.) m_vtx2Daug3MassHypo = gendata->particleMass(MC::PIPLUS).value();
56 if(m_vtx2Daug4MassHypo < 0.) m_vtx2Daug4MassHypo = gendata->particleMass(MC::PIPLUS).value();
57
58 return StatusCode::SUCCESS;
59 }
60
61 StatusCode PsiPlusPsiCascade::addBranches(const EventContext& ctx) const {
62 if ((m_vtx1Daug_num != 3 && m_vtx1Daug_num != 4) || (m_vtx2Daug_num != 3 && m_vtx2Daug_num != 4)) {
63 ATH_MSG_FATAL("Incorrect number of Psi daughters (should be 3 or 4)");
64 return StatusCode::FAILURE;
65 }
66
67 constexpr int topoN = 3;
68 if(m_cascadeOutputsKeys.size() != topoN) {
69 ATH_MSG_FATAL("Incorrect number of VtxContainers");
70 return StatusCode::FAILURE;
71 }
72 std::array<SG::WriteHandle<xAOD::VertexContainer>, topoN> VtxWriteHandles; int ikey(0);
74 VtxWriteHandles[ikey] = SG::WriteHandle<xAOD::VertexContainer>(key, ctx);
75 ATH_CHECK( VtxWriteHandles[ikey].record(std::make_unique<xAOD::VertexContainer>(), std::make_unique<xAOD::VertexAuxContainer>()) );
76 ikey++;
77 }
78
79 //----------------------------------------------------
80 // retrieve primary vertices
81 //----------------------------------------------------
83 ATH_CHECK( pvContainer.isValid() );
84 if (pvContainer.cptr()->size()==0) {
85 ATH_MSG_WARNING("You have no primary vertices: " << pvContainer.cptr()->size());
86 return StatusCode::RECOVERABLE;
87 }
88
89 //----------------------------------------------------
90 // Record refitted primary vertices
91 //----------------------------------------------------
93 if(m_refitPV) {
95 ATH_CHECK( refPvContainer.record(std::make_unique<xAOD::VertexContainer>(), std::make_unique<xAOD::VertexAuxContainer>()) );
96 }
97
98 std::vector<Trk::VxCascadeInfo*> cascadeinfoContainer;
99 std::vector<Trk::VxCascadeInfo*> cascadeinfoContainer_noConstr;
100 ATH_CHECK(performSearch(&cascadeinfoContainer,&cascadeinfoContainer_noConstr,ctx));
101
103 ATH_CHECK( evt.isValid() );
104 BPhysPVCascadeTools helper(&(*m_CascadeTools), evt.cptr());
105 helper.SetMinNTracksInPV(m_PV_minNTracks);
106
107 // Decorators for the main vertex: chi2, ndf, pt and pt error, plus the V0 vertex variables
108 SG::Decorator<VertexLinkVector> CascadeLinksDecor("CascadeVertexLinks");
109 SG::Decorator<VertexLinkVector> Psi1LinksDecor("Psi1VertexLinks");
110 SG::Decorator<VertexLinkVector> Psi2LinksDecor("Psi2VertexLinks");
111 SG::Decorator<float> chi2_decor("ChiSquared");
112 SG::Decorator<int> ndof_decor("nDoF");
113 SG::Decorator<float> chi2_nc_decor("ChiSquared_nc");
114 SG::Decorator<int> ndof_nc_decor("nDoF_nc");
115 SG::Decorator<float> Pt_decor("Pt");
116 SG::Decorator<float> PtErr_decor("PtErr");
117 SG::Decorator<float> chi2_SV1_decor("ChiSquared_SV1");
118 SG::Decorator<float> chi2_nc_SV1_decor("ChiSquared_nc_SV1");
119 SG::Decorator<float> chi2_V1_decor("ChiSquared_V1");
120 SG::Decorator<int> ndof_V1_decor("nDoF_V1");
121 SG::Decorator<float> lxy_SV1_decor("lxy_SV1");
122 SG::Decorator<float> lxyErr_SV1_decor("lxyErr_SV1");
123 SG::Decorator<float> a0xy_SV1_decor("a0xy_SV1");
124 SG::Decorator<float> a0xyErr_SV1_decor("a0xyErr_SV1");
125 SG::Decorator<float> a0z_SV1_decor("a0z_SV1");
126 SG::Decorator<float> a0zErr_SV1_decor("a0zErr_SV1");
127 SG::Decorator<float> chi2_SV2_decor("ChiSquared_SV2");
128 SG::Decorator<float> chi2_nc_SV2_decor("ChiSquared_nc_SV2");
129 SG::Decorator<float> chi2_V2_decor("ChiSquared_V2");
130 SG::Decorator<int> ndof_V2_decor("nDoF_V2");
131 SG::Decorator<float> lxy_SV2_decor("lxy_SV2");
132 SG::Decorator<float> lxyErr_SV2_decor("lxyErr_SV2");
133 SG::Decorator<float> a0xy_SV2_decor("a0xy_SV2");
134 SG::Decorator<float> a0xyErr_SV2_decor("a0xyErr_SV2");
135 SG::Decorator<float> a0z_SV2_decor("a0z_SV2");
136 SG::Decorator<float> a0zErr_SV2_decor("a0zErr_SV2");
137
138 // Get the container and identify the input Psi's
140 ATH_CHECK( psi1Container.isValid() );
142 ATH_CHECK( psi2Container.isValid() );
143
144 for(size_t ic=0; ic<cascadeinfoContainer.size(); ic++) {
145 Trk::VxCascadeInfo* cascade_info = cascadeinfoContainer[ic];
146 if(cascade_info==nullptr) {
147 ATH_MSG_ERROR("CascadeInfo is null");
148 continue;
149 }
150
151 Trk::VxCascadeInfo* cascade_info_noConstr = cascadeinfoContainer_noConstr[ic];
152
153 const std::vector<xAOD::Vertex*> &cascadeVertices = cascade_info->vertices();
154 if(cascadeVertices.size() != topoN) ATH_MSG_ERROR("Incorrect number of vertices");
155 if(cascadeVertices[0]==nullptr || cascadeVertices[1]==nullptr || cascadeVertices[2]==nullptr) ATH_MSG_ERROR("Error null vertex");
156 // Keep vertices
157 for(int i=0; i<topoN; i++) VtxWriteHandles[i].ptr()->push_back(cascadeVertices[i]);
158
159 cascade_info->setSVOwnership(false); // Prevent Container from deleting vertices
160 const auto mainVertex = cascadeVertices[2]; // this is the mother vertex
161 const std::vector< std::vector<TLorentzVector> > &moms = cascade_info->getParticleMoms();
162
163 // Set links to cascade vertices
164 std::vector<VertexLink> precedingVertexLinks;
165 VertexLink vertexLink1;
166 vertexLink1.setElement(cascadeVertices[0]);
167 vertexLink1.setStorableObject(*VtxWriteHandles[0].ptr());
168 if( vertexLink1.isValid() ) precedingVertexLinks.push_back( vertexLink1 );
169 VertexLink vertexLink2;
170 vertexLink2.setElement(cascadeVertices[1]);
171 vertexLink2.setStorableObject(*VtxWriteHandles[1].ptr());
172 if( vertexLink2.isValid() ) precedingVertexLinks.push_back( vertexLink2 );
173 CascadeLinksDecor(*mainVertex) = std::move(precedingVertexLinks);
174
175 // Identify the input Psi2
176 const xAOD::Vertex* psi2Vertex(0);
177 if(m_vtx2Daug_num==4) psi2Vertex = BPhysPVCascadeTools::FindVertex<4>(psi2Container.cptr(), cascadeVertices[1]);
178 else psi2Vertex = BPhysPVCascadeTools::FindVertex<3>(psi2Container.cptr(), cascadeVertices[1]);
179 // Identify the input Psi1
180 const xAOD::Vertex* psi1Vertex(0);
181 if(m_vtx1Daug_num==4) psi1Vertex = BPhysPVCascadeTools::FindVertex<4>(psi1Container.cptr(), cascadeVertices[0]);
182 else psi1Vertex = BPhysPVCascadeTools::FindVertex<3>(psi1Container.cptr(), cascadeVertices[0]);
183
184 // Set links to input vertices
185 std::vector<const xAOD::Vertex*> psi2VerticestoLink;
186 if(psi2Vertex) psi2VerticestoLink.push_back(psi2Vertex);
187 else ATH_MSG_WARNING("Could not find linking Jpsi");
188 if(!BPhysPVCascadeTools::LinkVertices(Psi2LinksDecor, psi2VerticestoLink, psi2Container.cptr(), mainVertex)) ATH_MSG_ERROR("Error decorating with Psi2 vertex");
189
190 std::vector<const xAOD::Vertex*> psi1VerticestoLink;
191 if(psi1Vertex) psi1VerticestoLink.push_back(psi1Vertex);
192 else ATH_MSG_WARNING("Could not find linking Psi1");
193 if(!BPhysPVCascadeTools::LinkVertices(Psi1LinksDecor, psi1VerticestoLink, psi1Container.cptr(), mainVertex)) ATH_MSG_ERROR("Error decorating with Psi1 vertex");
194
195 xAOD::BPhysHypoHelper vtx(m_hypoName, mainVertex);
196
197 // Get refitted track momenta from all vertices, charged tracks only
198 BPhysPVCascadeTools::SetVectorInfo(vtx, cascade_info);
199 vtx.setPass(true);
200
201 //
202 // Decorate main vertex
203 //
204 // mass, mass error
205 // https://gitlab.cern.ch/atlas/athena/-/blob/21.2/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/VxCascadeInfo.h
206 BPHYS_CHECK( vtx.setMass(m_CascadeTools->invariantMass(moms[2])) );
207 BPHYS_CHECK( vtx.setMassErr(m_CascadeTools->invariantMassError(moms[2],cascade_info->getCovariance()[1])) );
208 // pt and pT error (the default pt of mainVertex is != the pt of the full cascade fit!)
209 Pt_decor(*mainVertex) = m_CascadeTools->pT(moms[2]);
210 PtErr_decor(*mainVertex) = m_CascadeTools->pTError(moms[2],cascade_info->getCovariance()[1]);
211 // chi2 and ndof (the default chi2 of mainVertex is != the chi2 of the full cascade fit!)
212 chi2_decor(*mainVertex) = cascade_info->fitChi2();
213 ndof_decor(*mainVertex) = cascade_info->nDoF();
214 chi2_nc_decor(*mainVertex) = cascade_info_noConstr ? cascade_info_noConstr->fitChi2() : -999999.;
215 ndof_nc_decor(*mainVertex) = cascade_info_noConstr ? cascade_info_noConstr->nDoF() : -1;
216
217 // decorate the Psi1 vertex
218 chi2_SV1_decor(*cascadeVertices[0]) = m_V0Tools->chisq(cascadeVertices[0]);
219 chi2_nc_SV1_decor(*cascadeVertices[0]) = cascade_info_noConstr ? m_V0Tools->chisq(cascade_info_noConstr->vertices()[0]) : -999999.;
220 //
221 auto ndof = psi1Vertex ? m_V0Tools->ndof(psi1Vertex) : -999999.;
222 auto chisq = psi1Vertex ? m_V0Tools->chisq(psi1Vertex) : 999999.;
223 chi2_V1_decor(*cascadeVertices[0]) = chisq;
224 ndof_V1_decor(*cascadeVertices[0]) = ndof;
225 //
226 lxy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],mainVertex);
227 lxyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxyError(moms[0],cascade_info->getCovariance()[0],cascadeVertices[0],mainVertex);
228 a0z_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0z(moms[0],cascadeVertices[0],mainVertex);
229 a0zErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0zError(moms[0],cascade_info->getCovariance()[0],cascadeVertices[0],mainVertex);
230 a0xy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xy(moms[0],cascadeVertices[0],mainVertex);
231 a0xyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xyError(moms[0],cascade_info->getCovariance()[0],cascadeVertices[0],mainVertex);
232
233 // decorate the Psi2 vertex
234 chi2_SV2_decor(*cascadeVertices[1]) = m_V0Tools->chisq(cascadeVertices[1]);
235 chi2_nc_SV2_decor(*cascadeVertices[1]) = cascade_info_noConstr ? m_V0Tools->chisq(cascade_info_noConstr->vertices()[1]) : -999999.;
236 //
237 ndof = psi2Vertex ? m_V0Tools->ndof(psi1Vertex) : -999999.;
238 chisq = psi2Vertex ? m_V0Tools->chisq(psi1Vertex) : 999999.;
239 chi2_V2_decor(*cascadeVertices[1]) = chisq;
240 ndof_V2_decor(*cascadeVertices[1]) = ndof;
241 //
242 lxy_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->lxy(moms[1],cascadeVertices[1],mainVertex);
243 lxyErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->lxyError(moms[1],cascade_info->getCovariance()[1],cascadeVertices[1],mainVertex);
244 a0z_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0z(moms[1],cascadeVertices[1],mainVertex);
245 a0zErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0zError(moms[1],cascade_info->getCovariance()[1],cascadeVertices[1],mainVertex);
246 a0xy_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0xy(moms[1],cascadeVertices[1],mainVertex);
247 a0xyErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0xyError(moms[1],cascade_info->getCovariance()[1],cascadeVertices[1],mainVertex);
248
249 double Mass_Moth = m_CascadeTools->invariantMass(moms[2]); // size=2
250 ATH_CHECK(helper.FillCandwithRefittedVertices(m_refitPV, pvContainer.cptr(), m_refitPV ? refPvContainer.ptr() : 0, &(*m_pvRefitter), m_PV_max, m_DoVertexType, cascade_info, 2, Mass_Moth, vtx));
251 } // loop over cascadeinfoContainer
252
253 // Deleting cascadeinfo since this won't be stored.
254 // Vertices have been kept in m_cascadeOutputs and should be owned by their container
255 for (auto cascade_info : cascadeinfoContainer) delete cascade_info;
256 for (auto cascade_info_noConstr : cascadeinfoContainer_noConstr) delete cascade_info_noConstr;
257
258 return StatusCode::SUCCESS;
259 }
260
261 PsiPlusPsiCascade::PsiPlusPsiCascade(const std::string& type, const std::string& name, const IInterface* parent) : base_class(type,name,parent),
264 m_cascadeOutputsKeys({"PsiPlusPsiCascadeVtx1", "PsiPlusPsiCascadeVtx2", "PsiPlusPsiCascadeVtx3"}),
265 m_VxPrimaryCandidateName("PrimaryVertices"),
266 m_trackContainerName("InDetTrackParticles"),
267 m_eventInfo_key("EventInfo"),
268 m_jpsi1MassLower(0.0),
269 m_jpsi1MassUpper(20000.0),
270 m_jpsi2MassLower(0.0),
271 m_jpsi2MassUpper(20000.0),
272 m_diTrack1MassLower(-1.0),
273 m_diTrack1MassUpper(-1.0),
274 m_diTrack2MassLower(-1.0),
275 m_diTrack2MassUpper(-1.0),
276 m_psi1MassLower(0.0),
277 m_psi1MassUpper(25000.0),
278 m_psi2MassLower(0.0),
279 m_psi2MassUpper(25000.0),
280 m_MassLower(0.0),
281 m_MassUpper(31000.0),
282 m_vtx1Daug_num(4),
283 m_vtx1Daug1MassHypo(-1),
284 m_vtx1Daug2MassHypo(-1),
285 m_vtx1Daug3MassHypo(-1),
286 m_vtx1Daug4MassHypo(-1),
287 m_vtx2Daug_num(4),
288 m_vtx2Daug1MassHypo(-1),
289 m_vtx2Daug2MassHypo(-1),
290 m_vtx2Daug3MassHypo(-1),
291 m_vtx2Daug4MassHypo(-1),
292 m_massPsi1(-1),
293 m_massPsi2(-1),
294 m_massJpsi1(-1),
295 m_massJpsi2(-1),
296 m_massDiTrk1(-1),
297 m_massDiTrk2(-1),
298 m_constrPsi1(false),
299 m_constrPsi2(false),
300 m_constrJpsi1(false),
301 m_constrJpsi2(false),
302 m_constrDiTrk1(false),
303 m_constrDiTrk2(false),
304 m_chi2cut_Psi1(-1.0),
305 m_chi2cut_Psi2(-1.0),
306 m_chi2cut(-1.0),
307 m_removeDuplicatePairs(false),
308 m_maxCandidates(0),
309 m_iVertexFitter("Trk::TrkVKalVrtFitter"),
310 m_pvRefitter("Analysis::PrimaryVertexRefitter"),
311 m_V0Tools("Trk::V0Tools"),
312 m_CascadeTools("DerivationFramework::CascadeTools")
313 {
314 declareProperty("Psi1Vertices", m_vertexPsi1ContainerKey);
315 declareProperty("Psi2Vertices", m_vertexPsi2ContainerKey);
316 declareProperty("Psi1VtxHypoNames", m_vertexPsi1HypoNames);
317 declareProperty("Psi2VtxHypoNames", m_vertexPsi2HypoNames);
318 declareProperty("VxPrimaryCandidateName", m_VxPrimaryCandidateName);
319 declareProperty("TrackContainerName", m_trackContainerName);
320 declareProperty("RefPVContainerName", m_refPVContainerName = "RefittedPrimaryVertices");
321 declareProperty("Jpsi1MassLowerCut", m_jpsi1MassLower);
322 declareProperty("Jpsi1MassUpperCut", m_jpsi1MassUpper);
323 declareProperty("Jpsi2MassLowerCut", m_jpsi2MassLower);
324 declareProperty("Jpsi2MassUpperCut", m_jpsi2MassUpper);
325 declareProperty("DiTrack1MassLower", m_diTrack1MassLower); // only effective when m_vtx1Daug_num=4
326 declareProperty("DiTrack1MassUpper", m_diTrack1MassUpper); // only effective when m_vtx1Daug_num=4
327 declareProperty("DiTrack2MassLower", m_diTrack2MassLower); // only effective when m_vtx2Daug_num=4
328 declareProperty("DiTrack2MassUpper", m_diTrack2MassUpper); // only effective when m_vtx2Daug_num=4
329 declareProperty("Psi1MassLowerCut", m_psi1MassLower);
330 declareProperty("Psi1MassUpperCut", m_psi1MassUpper);
331 declareProperty("Psi2MassLowerCut", m_psi2MassLower);
332 declareProperty("Psi2MassUpperCut", m_psi2MassUpper);
333 declareProperty("MassLowerCut", m_MassLower);
334 declareProperty("MassUpperCut", m_MassUpper);
335 declareProperty("HypothesisName", m_hypoName = "TQ");
336 declareProperty("NumberOfPsi1Daughters", m_vtx1Daug_num); // 3 or 4 only
337 declareProperty("Vtx1Daug1MassHypo", m_vtx1Daug1MassHypo);
338 declareProperty("Vtx1Daug2MassHypo", m_vtx1Daug2MassHypo);
339 declareProperty("Vtx1Daug3MassHypo", m_vtx1Daug3MassHypo);
340 declareProperty("Vtx1Daug4MassHypo", m_vtx1Daug4MassHypo);
341 declareProperty("NumberOfPsi2Daughters", m_vtx2Daug_num); // 3 or 4 only
342 declareProperty("Vtx2Daug1MassHypo", m_vtx2Daug1MassHypo);
343 declareProperty("Vtx2Daug2MassHypo", m_vtx2Daug2MassHypo);
344 declareProperty("Vtx2Daug3MassHypo", m_vtx2Daug3MassHypo);
345 declareProperty("Vtx2Daug4MassHypo", m_vtx2Daug4MassHypo);
346 declareProperty("Jpsi1Mass", m_massJpsi1);
347 declareProperty("Jpsi2Mass", m_massJpsi2);
348 declareProperty("DiTrack1Mass", m_massDiTrk1);
349 declareProperty("DiTrack2Mass", m_massDiTrk2);
350 declareProperty("Psi1Mass", m_massPsi1);
351 declareProperty("Psi2Mass", m_massPsi2);
352 declareProperty("ApplyPsi1MassConstraint", m_constrPsi1);
353 declareProperty("ApplyPsi2MassConstraint", m_constrPsi2);
354 declareProperty("ApplyJpsi1MassConstraint", m_constrJpsi1);
355 declareProperty("ApplyJpsi2MassConstraint", m_constrJpsi2);
356 declareProperty("ApplyDiTrk1MassConstraint",m_constrDiTrk1); // only effective when m_vtx1Daug_num=4
357 declareProperty("ApplyDiTrk2MassConstraint",m_constrDiTrk2); // only effective when m_vtx2Daug_num=4
358 declareProperty("Chi2CutPsi1", m_chi2cut_Psi1);
359 declareProperty("Chi2CutPsi2", m_chi2cut_Psi2);
360 declareProperty("Chi2Cut", m_chi2cut);
361 declareProperty("RemoveDuplicatePairs", m_removeDuplicatePairs); // only effective when m_vertexPsi1ContainerKey == m_vertexPsi2ContainerKey
362 declareProperty("MaxCandidates", m_maxCandidates);
363 declareProperty("RefitPV", m_refitPV = true);
364 declareProperty("MaxnPV", m_PV_max = 1000);
365 declareProperty("MinNTracksInPV", m_PV_minNTracks = 0);
366 declareProperty("DoVertexType", m_DoVertexType = 7);
367 declareProperty("TrkVertexFitterTool", m_iVertexFitter);
368 declareProperty("PVRefitter", m_pvRefitter);
369 declareProperty("V0Tools", m_V0Tools);
370 declareProperty("CascadeTools", m_CascadeTools);
371 declareProperty("CascadeVertexCollections", m_cascadeOutputsKeys);
372 }
373
374 StatusCode PsiPlusPsiCascade::performSearch(std::vector<Trk::VxCascadeInfo*> *cascadeinfoContainer, std::vector<Trk::VxCascadeInfo*> *cascadeinfoContainer_noConstr, const EventContext& ctx) const {
375 ATH_MSG_DEBUG( "PsiPlusPsiCascade::performSearch" );
376 assert(cascadeinfoContainer!=nullptr && cascadeinfoContainer_noConstr!=nullptr);
377
378 // Get TrackParticle container (for setting links to the original tracks)
380 ATH_CHECK( trackContainer.isValid() );
381
382 std::vector<const xAOD::TrackParticle*> tracksJpsi1;
383 std::vector<const xAOD::TrackParticle*> tracksJpsi2;
384 std::vector<const xAOD::TrackParticle*> tracksDiTrk1;
385 std::vector<const xAOD::TrackParticle*> tracksDiTrk2;
386 std::vector<const xAOD::TrackParticle*> tracksPsi1;
387 std::vector<const xAOD::TrackParticle*> tracksPsi2;
388 std::vector<double> massesPsi1;
389 massesPsi1.push_back(m_vtx1Daug1MassHypo);
390 massesPsi1.push_back(m_vtx1Daug2MassHypo);
391 massesPsi1.push_back(m_vtx1Daug3MassHypo);
392 if(m_vtx1Daug_num==4) massesPsi1.push_back(m_vtx1Daug4MassHypo);
393 std::vector<double> massesPsi2;
394 massesPsi2.push_back(m_vtx2Daug1MassHypo);
395 massesPsi2.push_back(m_vtx2Daug2MassHypo);
396 massesPsi2.push_back(m_vtx2Daug3MassHypo);
397 if(m_vtx2Daug_num==4) massesPsi2.push_back(m_vtx2Daug4MassHypo);
398
399 // Get Psi1 container
401 ATH_CHECK( psi1Container.isValid() );
402
403 // Get Psi2 container
405 ATH_CHECK( psi2Container.isValid() );
406
407 // Select the Psi2 candidates before calling cascade fit
408 std::vector<const xAOD::Vertex*> selectedPsi2Candidates;
409 for(auto vxcItr=psi2Container.cptr()->cbegin(); vxcItr!=psi2Container.cptr()->cend(); ++vxcItr) {
410 // Check the passed flag first
411 const xAOD::Vertex* vtx = *vxcItr;
412 bool passed = false;
413 for(size_t i=0; i<m_vertexPsi2HypoNames.size(); i++) {
414 SG::Accessor<Char_t> flagAcc("passed_"+m_vertexPsi2HypoNames[i]);
415 if(flagAcc.isAvailable(*vtx) && flagAcc(*vtx)) {
416 passed |= 1;
417 }
418 }
419 if(m_vertexPsi2HypoNames.size() && !passed) continue;
420
421 // Check Psi2 candidate invariant mass and skip if need be
422 double mass_psi2 = m_V0Tools->invariantMass(*vxcItr, massesPsi2);
423 if (mass_psi2 < m_psi2MassLower || mass_psi2 > m_psi2MassUpper) continue;
424
425 // Add loose cut on Jpsi2 mass from Psi2 -> Jpsi2 pi+ pi-, or on phi mass from Ds+ -> phi pi+
426 TLorentzVector p4_mu1, p4_mu2;
427 p4_mu1.SetPtEtaPhiM( (*vxcItr)->trackParticle(0)->pt(),
428 (*vxcItr)->trackParticle(0)->eta(),
429 (*vxcItr)->trackParticle(0)->phi(), m_vtx2Daug1MassHypo);
430 p4_mu2.SetPtEtaPhiM( (*vxcItr)->trackParticle(1)->pt(),
431 (*vxcItr)->trackParticle(1)->eta(),
432 (*vxcItr)->trackParticle(1)->phi(), m_vtx2Daug2MassHypo);
433 double mass_jpsi2 = (p4_mu1 + p4_mu2).M();
434 if (mass_jpsi2 < m_jpsi2MassLower || mass_jpsi2 > m_jpsi2MassUpper) continue;
435
437 TLorentzVector p4_trk1, p4_trk2;
438 p4_trk1.SetPtEtaPhiM( (*vxcItr)->trackParticle(2)->pt(),
439 (*vxcItr)->trackParticle(2)->eta(),
440 (*vxcItr)->trackParticle(2)->phi(), m_vtx2Daug3MassHypo);
441 p4_trk2.SetPtEtaPhiM( (*vxcItr)->trackParticle(3)->pt(),
442 (*vxcItr)->trackParticle(3)->eta(),
443 (*vxcItr)->trackParticle(3)->phi(), m_vtx2Daug4MassHypo);
444 double mass_diTrk2 = (p4_trk1 + p4_trk2).M();
445 if (mass_diTrk2 < m_diTrack2MassLower || mass_diTrk2 > m_diTrack2MassUpper) continue;
446 }
447
448 double chi2DOF = (*vxcItr)->chiSquared()/(*vxcItr)->numberDoF();
449 if(m_chi2cut_Psi2>0 && chi2DOF>m_chi2cut_Psi2) continue;
450
451 selectedPsi2Candidates.push_back(*vxcItr);
452 }
453 if(selectedPsi2Candidates.size()==0) return StatusCode::SUCCESS;
454
455 // Select the Psi1 candidates before calling cascade fit
456 std::vector<const xAOD::Vertex*> selectedPsi1Candidates;
457 for(auto vxcItr=psi1Container.cptr()->cbegin(); vxcItr!=psi1Container.cptr()->cend(); ++vxcItr) {
458 // Check the passed flag first
459 const xAOD::Vertex* vtx = *vxcItr;
460 bool passed = false;
461 for(size_t i=0; i<m_vertexPsi1HypoNames.size(); i++) {
462 SG::Accessor<Char_t> flagAcc("passed_"+m_vertexPsi1HypoNames[i]);
463 if(flagAcc.isAvailable(*vtx) && flagAcc(*vtx)) {
464 passed |= 1;
465 }
466 }
467 if(m_vertexPsi1HypoNames.size() && !passed) continue;
468
469 // Check Psi candidate invariant mass and skip if need be
470 double mass_psi1 = m_V0Tools->invariantMass(*vxcItr,massesPsi1);
471 if(mass_psi1 < m_psi1MassLower || mass_psi1 > m_psi1MassUpper) continue;
472
473 // Add loose cut on Jpsi1 mass from Psi1 -> Jpsi1 pi+ pi-, or on phi mass from Ds+ -> phi pi+
474 TLorentzVector p4_mu1, p4_mu2;
475 p4_mu1.SetPtEtaPhiM( (*vxcItr)->trackParticle(0)->pt(),
476 (*vxcItr)->trackParticle(0)->eta(),
477 (*vxcItr)->trackParticle(0)->phi(), m_vtx1Daug1MassHypo);
478 p4_mu2.SetPtEtaPhiM( (*vxcItr)->trackParticle(1)->pt(),
479 (*vxcItr)->trackParticle(1)->eta(),
480 (*vxcItr)->trackParticle(1)->phi(), m_vtx1Daug2MassHypo);
481 double mass_jpsi1 = (p4_mu1 + p4_mu2).M();
482 if (mass_jpsi1 < m_jpsi1MassLower || mass_jpsi1 > m_jpsi1MassUpper) continue;
483
485 TLorentzVector p4_trk1, p4_trk2;
486 p4_trk1.SetPtEtaPhiM( (*vxcItr)->trackParticle(2)->pt(),
487 (*vxcItr)->trackParticle(2)->eta(),
488 (*vxcItr)->trackParticle(2)->phi(), m_vtx1Daug3MassHypo);
489 p4_trk2.SetPtEtaPhiM( (*vxcItr)->trackParticle(3)->pt(),
490 (*vxcItr)->trackParticle(3)->eta(),
491 (*vxcItr)->trackParticle(3)->phi(), m_vtx1Daug4MassHypo);
492 double mass_diTrk1 = (p4_trk1 + p4_trk2).M();
493 if (mass_diTrk1 < m_diTrack1MassLower || mass_diTrk1 > m_diTrack1MassUpper) continue;
494 }
495
496 double chi2DOF = (*vxcItr)->chiSquared()/(*vxcItr)->numberDoF();
497 if(m_chi2cut_Psi1>0 && chi2DOF>m_chi2cut_Psi1) continue;
498
499 selectedPsi1Candidates.push_back(*vxcItr);
500 }
501 if(selectedPsi1Candidates.size()==0) return StatusCode::SUCCESS;
502
503 std::vector<std::pair<const xAOD::Vertex*, const xAOD::Vertex*> > candidatePairs;
504 for(auto psi1Itr=selectedPsi1Candidates.cbegin(); psi1Itr!=selectedPsi1Candidates.cend(); ++psi1Itr) {
505 tracksPsi1.clear();
506 for(size_t i=0; i<(*psi1Itr)->nTrackParticles(); i++) tracksPsi1.push_back((*psi1Itr)->trackParticle(i));
507 for(auto psi2Itr=selectedPsi2Candidates.cbegin(); psi2Itr!=selectedPsi2Candidates.cend(); ++psi2Itr) {
508 bool skip = false;
509 for(size_t j=0; j<(*psi2Itr)->nTrackParticles(); j++) {
510 if(std::find(tracksPsi1.cbegin(), tracksPsi1.cend(), (*psi2Itr)->trackParticle(j)) != tracksPsi1.cend()) { skip = true; break; }
511 }
512 if(skip) continue;
514 for(size_t ic=0; ic<candidatePairs.size(); ic++) {
515 const xAOD::Vertex* psi1Vertex = candidatePairs[ic].first;
516 const xAOD::Vertex* psi2Vertex = candidatePairs[ic].second;
517 if((psi1Vertex == *psi1Itr && psi2Vertex == *psi2Itr) || (psi1Vertex == *psi2Itr && psi2Vertex == *psi1Itr)) { skip = true; break; }
518 }
519 }
520 if(skip) continue;
521 candidatePairs.push_back(std::pair<const xAOD::Vertex*, const xAOD::Vertex*>(*psi1Itr,*psi2Itr));
522 }
523 }
524
525 std::sort( candidatePairs.begin(), candidatePairs.end(), [](std::pair<const xAOD::Vertex*, const xAOD::Vertex*> a, std::pair<const xAOD::Vertex*, const xAOD::Vertex*> b) { return a.first->chiSquared()/a.first->numberDoF()+a.second->chiSquared()/a.second->numberDoF() < b.first->chiSquared()/b.first->numberDoF()+b.second->chiSquared()/b.second->numberDoF(); } );
526 if(m_maxCandidates>0 && candidatePairs.size()>m_maxCandidates) {
527 candidatePairs.erase(candidatePairs.begin()+m_maxCandidates, candidatePairs.end());
528 }
529
530 for(size_t ic=0; ic<candidatePairs.size(); ic++) {
531 const xAOD::Vertex* psi1Vertex = candidatePairs[ic].first;
532 const xAOD::Vertex* psi2Vertex = candidatePairs[ic].second;
533
534 tracksPsi1.clear();
535 for(size_t it=0; it<psi1Vertex->nTrackParticles(); it++) tracksPsi1.push_back(psi1Vertex->trackParticle(it));
536 if (tracksPsi1.size() != massesPsi1.size()) {
537 ATH_MSG_ERROR("Problems with Psi1 input: number of tracks or track mass inputs is not correct!");
538 }
539 tracksPsi2.clear();
540 for(size_t it=0; it<psi2Vertex->nTrackParticles(); it++) tracksPsi2.push_back(psi2Vertex->trackParticle(it));
541 if (tracksPsi2.size() != massesPsi2.size()) {
542 ATH_MSG_ERROR("Problems with Psi2 input: number of tracks or track mass inputs is not correct!");
543 }
544
545 tracksJpsi1.clear();
546 tracksJpsi1.push_back(psi1Vertex->trackParticle(0));
547 tracksJpsi1.push_back(psi1Vertex->trackParticle(1));
548 tracksDiTrk1.clear();
549 if(m_vtx1Daug_num==4) {
550 tracksDiTrk1.push_back(psi1Vertex->trackParticle(2));
551 tracksDiTrk1.push_back(psi1Vertex->trackParticle(3));
552 }
553 tracksJpsi2.clear();
554 tracksJpsi2.push_back(psi2Vertex->trackParticle(0));
555 tracksJpsi2.push_back(psi2Vertex->trackParticle(1));
556 tracksDiTrk2.clear();
557 if(m_vtx2Daug_num==4) {
558 tracksDiTrk2.push_back(psi2Vertex->trackParticle(2));
559 tracksDiTrk2.push_back(psi2Vertex->trackParticle(3));
560 }
561
562 TLorentzVector p4_moth;
563 TLorentzVector tmp;
564 for(size_t it=0; it<psi1Vertex->nTrackParticles(); it++) {
565 tmp.SetPtEtaPhiM(psi1Vertex->trackParticle(it)->pt(),psi1Vertex->trackParticle(it)->eta(),psi1Vertex->trackParticle(it)->phi(),massesPsi1[it]);
566 p4_moth += tmp;
567 }
568 for(size_t it=0; it<psi2Vertex->nTrackParticles(); it++) {
569 tmp.SetPtEtaPhiM(psi2Vertex->trackParticle(it)->pt(),psi2Vertex->trackParticle(it)->eta(),psi2Vertex->trackParticle(it)->phi(),massesPsi2[it]);
570 p4_moth += tmp;
571 }
572 if (p4_moth.M() < m_MassLower || p4_moth.M() > m_MassUpper) continue;
573
574 // Apply the user's settings to the fitter
575 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
576 // Robustness: http://cdsweb.cern.ch/record/685551
577 int robustness = 0;
578 m_iVertexFitter->setRobustness(robustness, *state);
579 // Build up the topology
580 // Vertex list
581 std::vector<Trk::VertexID> vrtList;
582 // Psi1 vertex
583 Trk::VertexID vID1;
584 // https://gitlab.cern.ch/atlas/athena/-/blob/21.2/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
585 if (m_constrPsi1) {
586 vID1 = m_iVertexFitter->startVertex(tracksPsi1,massesPsi1,*state,m_massPsi1);
587 } else {
588 vID1 = m_iVertexFitter->startVertex(tracksPsi1,massesPsi1,*state);
589 }
590 vrtList.push_back(vID1);
591 // Psi2 vertex
592 Trk::VertexID vID2;
593 if (m_constrPsi2) {
594 vID2 = m_iVertexFitter->nextVertex(tracksPsi2,massesPsi2,*state,m_massPsi2);
595 } else {
596 vID2 = m_iVertexFitter->nextVertex(tracksPsi2,massesPsi2,*state);
597 }
598 vrtList.push_back(vID2);
599 // Mother vertex including Psi1 and Psi2
600 std::vector<const xAOD::TrackParticle*> tp; tp.clear();
601 std::vector<double> tp_masses; tp_masses.clear();
602 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList,*state);
603 if (m_constrJpsi1) {
604 std::vector<Trk::VertexID> cnstV; cnstV.clear();
605 if ( !m_iVertexFitter->addMassConstraint(vID1,tracksJpsi1,cnstV,*state,m_massJpsi1).isSuccess() ) {
606 ATH_MSG_WARNING("addMassConstraint for Jpsi1 failed");
607 }
608 }
610 std::vector<Trk::VertexID> cnstV; cnstV.clear();
611 if ( !m_iVertexFitter->addMassConstraint(vID1,tracksDiTrk1,cnstV,*state,m_massDiTrk1).isSuccess() ) {
612 ATH_MSG_WARNING("addMassConstraint for DiTrk1 failed");
613 }
614 }
615 if (m_constrJpsi2) {
616 std::vector<Trk::VertexID> cnstV; cnstV.clear();
617 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJpsi2,cnstV,*state,m_massJpsi2).isSuccess() ) {
618 ATH_MSG_WARNING("addMassConstraint for Jpsi2 failed");
619 }
620 }
622 std::vector<Trk::VertexID> cnstV; cnstV.clear();
623 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksDiTrk2,cnstV,*state,m_massDiTrk2).isSuccess() ) {
624 ATH_MSG_WARNING("addMassConstraint for DiTrk2 failed");
625 }
626 }
627 // Do the work
628 std::unique_ptr<Trk::VxCascadeInfo> result(m_iVertexFitter->fitCascade(*state));
629
630 bool pass = false;
631 if (result != nullptr) {
632 for(auto v : result->vertices()) {
633 if(v->nTrackParticles()==0) {
634 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
635 v->setTrackParticleLinks(nullLinkVector);
636 }
637 }
638 // reset links to original tracks
639 BPhysPVCascadeTools::PrepareVertexLinks(result.get(), trackContainer.cptr());
640
641 // necessary to prevent memory leak
642 result->setSVOwnership(true);
643
644 // Chi2/DOF cut
645 double chi2DOF = result->fitChi2()/result->nDoF();
646 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
647
648 if(chi2CutPassed) {
649 cascadeinfoContainer->push_back(result.release());
650 pass = true;
651 }
652 }
653
654 // do cascade fit again without any mass constraints
655 if(pass) {
657 std::unique_ptr<Trk::IVKalState> state (m_iVertexFitter->makeState(ctx));
658 m_iVertexFitter->setRobustness(robustness, *state);
659 std::vector<Trk::VertexID> vrtList_nc;
660 // Psi1 vertex
661 Trk::VertexID vID1_nc = m_iVertexFitter->startVertex(tracksPsi1,massesPsi1,*state);
662 vrtList_nc.push_back(vID1_nc);
663 // Psi2 vertex
664 Trk::VertexID vID2_nc = m_iVertexFitter->nextVertex(tracksPsi2,massesPsi2,*state);
665 vrtList_nc.push_back(vID2_nc);
666 // Mother vertex including Psi1 and Psi2
667 std::vector<const xAOD::TrackParticle*> tp; tp.clear();
668 std::vector<double> tp_masses; tp_masses.clear();
669 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList_nc,*state);
670 // Do the work
671 std::unique_ptr<Trk::VxCascadeInfo> result_nc(m_iVertexFitter->fitCascade(*state));
672
673 if (result_nc != nullptr) {
674 for(auto v : result_nc->vertices()) {
675 if(v->nTrackParticles()==0) {
676 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
677 v->setTrackParticleLinks(nullLinkVector);
678 }
679 }
680 // reset links to original tracks
681 BPhysPVCascadeTools::PrepareVertexLinks(result_nc.get(), trackContainer.cptr());
682
683 // necessary to prevent memory leak
684 result_nc->setSVOwnership(true);
685 cascadeinfoContainer_noConstr->push_back(result_nc.release());
686 }
687 else cascadeinfoContainer_noConstr->push_back(0);
688 }
689 else cascadeinfoContainer_noConstr->push_back(0);
690 }
691 } //Iterate over candidatePairs
692
693 return StatusCode::SUCCESS;
694 }
695}
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_DEBUG(x,...)
#define ATH_MSG_ERROR(x,...)
#define ATH_MSG_WARNING(x,...)
#define ATH_MSG_FATAL(x,...)
#define BPHYS_CHECK(EXP)
Useful CHECK macro.
: B-physics xAOD helpers.
ATLAS-specific HepMC functions.
bool passed(DecisionID id, const DecisionIDContainer &)
checks if required decision ID is in the set of IDs in the container
static Double_t a
static void PrepareVertexLinks(Trk::VxCascadeInfo *result, const xAOD::TrackParticleContainer *importedTrackCollection)
static const xAOD::Vertex * FindVertex(const xAOD::VertexContainer *c, const xAOD::Vertex *v)
static void SetVectorInfo(xAOD::BPhysHelper &, const Trk::VxCascadeInfo *)
static bool LinkVertices(SG::Decorator< VertexLinkVector > &decor, const std::vector< const xAOD::Vertex * > &vertices, const xAOD::VertexContainer *vertexContainer, const xAOD::Vertex *vert)
SG::WriteHandleKeyArray< xAOD::VertexContainer > m_cascadeOutputsKeys
SG::ReadHandleKey< xAOD::VertexContainer > m_vertexPsi1ContainerKey
virtual StatusCode initialize() override
PublicToolHandle< DerivationFramework::CascadeTools > m_CascadeTools
SG::ReadHandleKey< xAOD::EventInfo > m_eventInfo_key
PublicToolHandle< Trk::V0Tools > m_V0Tools
ToolHandle< Trk::TrkVKalVrtFitter > m_iVertexFitter
std::vector< std::string > m_vertexPsi1HypoNames
virtual StatusCode addBranches(const EventContext &ctx) const override
PublicToolHandle< Analysis::PrimaryVertexRefitter > m_pvRefitter
SG::ReadHandleKey< xAOD::VertexContainer > m_vertexPsi2ContainerKey
std::vector< std::string > m_vertexPsi2HypoNames
StatusCode performSearch(std::vector< Trk::VxCascadeInfo * > *cascadeinfoContainer, std::vector< Trk::VxCascadeInfo * > *cascadeinfoContainer_noConstr, const EventContext &ctx) const
SG::ReadHandleKey< xAOD::VertexContainer > m_VxPrimaryCandidateName
Name of primary vertex container.
PsiPlusPsiCascade(const std::string &type, const std::string &name, const IInterface *parent)
SG::WriteHandleKey< xAOD::VertexContainer > m_refPVContainerName
SG::ReadHandleKey< xAOD::TrackParticleContainer > m_trackContainerName
Helper class to provide type-safe access to aux data.
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
virtual bool isValid() override final
Can the handle be successfully dereferenced?
const_pointer_type cptr()
Dereference the pointer.
Property holding a SG store/key/clid from which a WriteHandle is made.
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
pointer_type ptr()
Dereference the pointer.
double fitChi2() const
const std::vector< Amg::MatrixX > & getCovariance() const
const std::vector< std::vector< TLorentzVector > > & getParticleMoms() const
const std::vector< xAOD::Vertex * > & vertices() const
void setSVOwnership(bool Ownership)
bool setMass(const float val)
Set given invariant mass and its error.
bool setPass(bool passVal)
get the pass flag for this hypothesis
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.
size_t nTrackParticles() const
Get the number of tracks associated with this vertex.
const TrackParticle * trackParticle(size_t i) const
Get the pointer to a given track that was used in vertex reco.
THE reconstruction tool.
ElementLink< xAOD::VertexContainer > VertexLink
std::vector< VertexLink > VertexLinkVector
static const int PSI2S
static const int MUON
static const int PIPLUS
static const int JPSI
SG::Decorator< T, ALLOC > Decorator
Helper class to provide type-safe access to aux data, specialized for JaggedVecElt.
Definition AuxElement.h:576
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
Vertex_v1 Vertex
Define the latest version of the vertex class.