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::AuxElement::Decorator<VertexLinkVector> CascadeLinksDecor("CascadeVertexLinks");
109 SG::AuxElement::Decorator<VertexLinkVector> Psi1LinksDecor("Psi1VertexLinks");
110 SG::AuxElement::Decorator<VertexLinkVector> Psi2LinksDecor("Psi2VertexLinks");
111 SG::AuxElement::Decorator<float> chi2_decor("ChiSquared");
112 SG::AuxElement::Decorator<int> ndof_decor("nDoF");
113 SG::AuxElement::Decorator<float> chi2_nc_decor("ChiSquared_nc");
114 SG::AuxElement::Decorator<int> ndof_nc_decor("nDoF_nc");
115 SG::AuxElement::Decorator<float> Pt_decor("Pt");
116 SG::AuxElement::Decorator<float> PtErr_decor("PtErr");
117 SG::AuxElement::Decorator<float> chi2_SV1_decor("ChiSquared_SV1");
118 SG::AuxElement::Decorator<float> chi2_nc_SV1_decor("ChiSquared_nc_SV1");
119 SG::AuxElement::Decorator<float> chi2_V1_decor("ChiSquared_V1");
120 SG::AuxElement::Decorator<int> ndof_V1_decor("nDoF_V1");
121 SG::AuxElement::Decorator<float> lxy_SV1_decor("lxy_SV1");
122 SG::AuxElement::Decorator<float> lxyErr_SV1_decor("lxyErr_SV1");
123 SG::AuxElement::Decorator<float> a0xy_SV1_decor("a0xy_SV1");
124 SG::AuxElement::Decorator<float> a0xyErr_SV1_decor("a0xyErr_SV1");
125 SG::AuxElement::Decorator<float> a0z_SV1_decor("a0z_SV1");
126 SG::AuxElement::Decorator<float> a0zErr_SV1_decor("a0zErr_SV1");
127 SG::AuxElement::Decorator<float> chi2_SV2_decor("ChiSquared_SV2");
128 SG::AuxElement::Decorator<float> chi2_nc_SV2_decor("ChiSquared_nc_SV2");
129 SG::AuxElement::Decorator<float> chi2_V2_decor("ChiSquared_V2");
130 SG::AuxElement::Decorator<int> ndof_V2_decor("nDoF_V2");
131 SG::AuxElement::Decorator<float> lxy_SV2_decor("lxy_SV2");
132 SG::AuxElement::Decorator<float> lxyErr_SV2_decor("lxyErr_SV2");
133 SG::AuxElement::Decorator<float> a0xy_SV2_decor("a0xy_SV2");
134 SG::AuxElement::Decorator<float> a0xyErr_SV2_decor("a0xyErr_SV2");
135 SG::AuxElement::Decorator<float> a0z_SV2_decor("a0z_SV2");
136 SG::AuxElement::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) = 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 chi2_V1_decor(*cascadeVertices[0]) = m_V0Tools->chisq(psi1Vertex);
221 ndof_V1_decor(*cascadeVertices[0]) = m_V0Tools->ndof(psi1Vertex);
222 lxy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],mainVertex);
223 lxyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxyError(moms[0],cascade_info->getCovariance()[0],cascadeVertices[0],mainVertex);
224 a0z_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0z(moms[0],cascadeVertices[0],mainVertex);
225 a0zErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0zError(moms[0],cascade_info->getCovariance()[0],cascadeVertices[0],mainVertex);
226 a0xy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xy(moms[0],cascadeVertices[0],mainVertex);
227 a0xyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xyError(moms[0],cascade_info->getCovariance()[0],cascadeVertices[0],mainVertex);
228
229 // decorate the Psi2 vertex
230 chi2_SV2_decor(*cascadeVertices[1]) = m_V0Tools->chisq(cascadeVertices[1]);
231 chi2_nc_SV2_decor(*cascadeVertices[1]) = cascade_info_noConstr ? m_V0Tools->chisq(cascade_info_noConstr->vertices()[1]) : -999999.;
232 chi2_V2_decor(*cascadeVertices[1]) = m_V0Tools->chisq(psi2Vertex);
233 ndof_V2_decor(*cascadeVertices[1]) = m_V0Tools->ndof(psi2Vertex);
234 lxy_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->lxy(moms[1],cascadeVertices[1],mainVertex);
235 lxyErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->lxyError(moms[1],cascade_info->getCovariance()[1],cascadeVertices[1],mainVertex);
236 a0z_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0z(moms[1],cascadeVertices[1],mainVertex);
237 a0zErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0zError(moms[1],cascade_info->getCovariance()[1],cascadeVertices[1],mainVertex);
238 a0xy_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0xy(moms[1],cascadeVertices[1],mainVertex);
239 a0xyErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0xyError(moms[1],cascade_info->getCovariance()[1],cascadeVertices[1],mainVertex);
240
241 double Mass_Moth = m_CascadeTools->invariantMass(moms[2]); // size=2
242 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));
243 } // loop over cascadeinfoContainer
244
245 // Deleting cascadeinfo since this won't be stored.
246 // Vertices have been kept in m_cascadeOutputs and should be owned by their container
247 for (auto cascade_info : cascadeinfoContainer) delete cascade_info;
248 for (auto cascade_info_noConstr : cascadeinfoContainer_noConstr) delete cascade_info_noConstr;
249
250 return StatusCode::SUCCESS;
251 }
252
253 PsiPlusPsiCascade::PsiPlusPsiCascade(const std::string& type, const std::string& name, const IInterface* parent) : base_class(type,name,parent),
256 m_cascadeOutputsKeys({"PsiPlusPsiCascadeVtx1", "PsiPlusPsiCascadeVtx2", "PsiPlusPsiCascadeVtx3"}),
257 m_VxPrimaryCandidateName("PrimaryVertices"),
258 m_trackContainerName("InDetTrackParticles"),
259 m_eventInfo_key("EventInfo"),
260 m_jpsi1MassLower(0.0),
261 m_jpsi1MassUpper(20000.0),
262 m_jpsi2MassLower(0.0),
263 m_jpsi2MassUpper(20000.0),
264 m_diTrack1MassLower(-1.0),
265 m_diTrack1MassUpper(-1.0),
266 m_diTrack2MassLower(-1.0),
267 m_diTrack2MassUpper(-1.0),
268 m_psi1MassLower(0.0),
269 m_psi1MassUpper(25000.0),
270 m_psi2MassLower(0.0),
271 m_psi2MassUpper(25000.0),
272 m_MassLower(0.0),
273 m_MassUpper(31000.0),
274 m_vtx1Daug_num(4),
275 m_vtx1Daug1MassHypo(-1),
276 m_vtx1Daug2MassHypo(-1),
277 m_vtx1Daug3MassHypo(-1),
278 m_vtx1Daug4MassHypo(-1),
279 m_vtx2Daug_num(4),
280 m_vtx2Daug1MassHypo(-1),
281 m_vtx2Daug2MassHypo(-1),
282 m_vtx2Daug3MassHypo(-1),
283 m_vtx2Daug4MassHypo(-1),
284 m_massPsi1(-1),
285 m_massPsi2(-1),
286 m_massJpsi1(-1),
287 m_massJpsi2(-1),
288 m_massDiTrk1(-1),
289 m_massDiTrk2(-1),
290 m_constrPsi1(false),
291 m_constrPsi2(false),
292 m_constrJpsi1(false),
293 m_constrJpsi2(false),
294 m_constrDiTrk1(false),
295 m_constrDiTrk2(false),
296 m_chi2cut_Psi1(-1.0),
297 m_chi2cut_Psi2(-1.0),
298 m_chi2cut(-1.0),
299 m_removeDuplicatePairs(false),
300 m_maxCandidates(0),
301 m_iVertexFitter("Trk::TrkVKalVrtFitter"),
302 m_pvRefitter("Analysis::PrimaryVertexRefitter"),
303 m_V0Tools("Trk::V0Tools"),
304 m_CascadeTools("DerivationFramework::CascadeTools")
305 {
306 declareProperty("Psi1Vertices", m_vertexPsi1ContainerKey);
307 declareProperty("Psi2Vertices", m_vertexPsi2ContainerKey);
308 declareProperty("Psi1VtxHypoNames", m_vertexPsi1HypoNames);
309 declareProperty("Psi2VtxHypoNames", m_vertexPsi2HypoNames);
310 declareProperty("VxPrimaryCandidateName", m_VxPrimaryCandidateName);
311 declareProperty("TrackContainerName", m_trackContainerName);
312 declareProperty("RefPVContainerName", m_refPVContainerName = "RefittedPrimaryVertices");
313 declareProperty("Jpsi1MassLowerCut", m_jpsi1MassLower);
314 declareProperty("Jpsi1MassUpperCut", m_jpsi1MassUpper);
315 declareProperty("Jpsi2MassLowerCut", m_jpsi2MassLower);
316 declareProperty("Jpsi2MassUpperCut", m_jpsi2MassUpper);
317 declareProperty("DiTrack1MassLower", m_diTrack1MassLower); // only effective when m_vtx1Daug_num=4
318 declareProperty("DiTrack1MassUpper", m_diTrack1MassUpper); // only effective when m_vtx1Daug_num=4
319 declareProperty("DiTrack2MassLower", m_diTrack2MassLower); // only effective when m_vtx2Daug_num=4
320 declareProperty("DiTrack2MassUpper", m_diTrack2MassUpper); // only effective when m_vtx2Daug_num=4
321 declareProperty("Psi1MassLowerCut", m_psi1MassLower);
322 declareProperty("Psi1MassUpperCut", m_psi1MassUpper);
323 declareProperty("Psi2MassLowerCut", m_psi2MassLower);
324 declareProperty("Psi2MassUpperCut", m_psi2MassUpper);
325 declareProperty("MassLowerCut", m_MassLower);
326 declareProperty("MassUpperCut", m_MassUpper);
327 declareProperty("HypothesisName", m_hypoName = "TQ");
328 declareProperty("NumberOfPsi1Daughters", m_vtx1Daug_num); // 3 or 4 only
329 declareProperty("Vtx1Daug1MassHypo", m_vtx1Daug1MassHypo);
330 declareProperty("Vtx1Daug2MassHypo", m_vtx1Daug2MassHypo);
331 declareProperty("Vtx1Daug3MassHypo", m_vtx1Daug3MassHypo);
332 declareProperty("Vtx1Daug4MassHypo", m_vtx1Daug4MassHypo);
333 declareProperty("NumberOfPsi2Daughters", m_vtx2Daug_num); // 3 or 4 only
334 declareProperty("Vtx2Daug1MassHypo", m_vtx2Daug1MassHypo);
335 declareProperty("Vtx2Daug2MassHypo", m_vtx2Daug2MassHypo);
336 declareProperty("Vtx2Daug3MassHypo", m_vtx2Daug3MassHypo);
337 declareProperty("Vtx2Daug4MassHypo", m_vtx2Daug4MassHypo);
338 declareProperty("Jpsi1Mass", m_massJpsi1);
339 declareProperty("Jpsi2Mass", m_massJpsi2);
340 declareProperty("DiTrack1Mass", m_massDiTrk1);
341 declareProperty("DiTrack2Mass", m_massDiTrk2);
342 declareProperty("Psi1Mass", m_massPsi1);
343 declareProperty("Psi2Mass", m_massPsi2);
344 declareProperty("ApplyPsi1MassConstraint", m_constrPsi1);
345 declareProperty("ApplyPsi2MassConstraint", m_constrPsi2);
346 declareProperty("ApplyJpsi1MassConstraint", m_constrJpsi1);
347 declareProperty("ApplyJpsi2MassConstraint", m_constrJpsi2);
348 declareProperty("ApplyDiTrk1MassConstraint",m_constrDiTrk1); // only effective when m_vtx1Daug_num=4
349 declareProperty("ApplyDiTrk2MassConstraint",m_constrDiTrk2); // only effective when m_vtx2Daug_num=4
350 declareProperty("Chi2CutPsi1", m_chi2cut_Psi1);
351 declareProperty("Chi2CutPsi2", m_chi2cut_Psi2);
352 declareProperty("Chi2Cut", m_chi2cut);
353 declareProperty("RemoveDuplicatePairs", m_removeDuplicatePairs); // only effective when m_vertexPsi1ContainerKey == m_vertexPsi2ContainerKey
354 declareProperty("MaxCandidates", m_maxCandidates);
355 declareProperty("RefitPV", m_refitPV = true);
356 declareProperty("MaxnPV", m_PV_max = 1000);
357 declareProperty("MinNTracksInPV", m_PV_minNTracks = 0);
358 declareProperty("DoVertexType", m_DoVertexType = 7);
359 declareProperty("TrkVertexFitterTool", m_iVertexFitter);
360 declareProperty("PVRefitter", m_pvRefitter);
361 declareProperty("V0Tools", m_V0Tools);
362 declareProperty("CascadeTools", m_CascadeTools);
363 declareProperty("CascadeVertexCollections", m_cascadeOutputsKeys);
364 }
365
366 StatusCode PsiPlusPsiCascade::performSearch(std::vector<Trk::VxCascadeInfo*> *cascadeinfoContainer, std::vector<Trk::VxCascadeInfo*> *cascadeinfoContainer_noConstr, const EventContext& ctx) const {
367 ATH_MSG_DEBUG( "PsiPlusPsiCascade::performSearch" );
368 assert(cascadeinfoContainer!=nullptr && cascadeinfoContainer_noConstr!=nullptr);
369
370 // Get TrackParticle container (for setting links to the original tracks)
372 ATH_CHECK( trackContainer.isValid() );
373
374 std::vector<const xAOD::TrackParticle*> tracksJpsi1;
375 std::vector<const xAOD::TrackParticle*> tracksJpsi2;
376 std::vector<const xAOD::TrackParticle*> tracksDiTrk1;
377 std::vector<const xAOD::TrackParticle*> tracksDiTrk2;
378 std::vector<const xAOD::TrackParticle*> tracksPsi1;
379 std::vector<const xAOD::TrackParticle*> tracksPsi2;
380 std::vector<double> massesPsi1;
381 massesPsi1.push_back(m_vtx1Daug1MassHypo);
382 massesPsi1.push_back(m_vtx1Daug2MassHypo);
383 massesPsi1.push_back(m_vtx1Daug3MassHypo);
384 if(m_vtx1Daug_num==4) massesPsi1.push_back(m_vtx1Daug4MassHypo);
385 std::vector<double> massesPsi2;
386 massesPsi2.push_back(m_vtx2Daug1MassHypo);
387 massesPsi2.push_back(m_vtx2Daug2MassHypo);
388 massesPsi2.push_back(m_vtx2Daug3MassHypo);
389 if(m_vtx2Daug_num==4) massesPsi2.push_back(m_vtx2Daug4MassHypo);
390
391 // Get Psi1 container
393 ATH_CHECK( psi1Container.isValid() );
394
395 // Get Psi2 container
397 ATH_CHECK( psi2Container.isValid() );
398
399 // Select the Psi2 candidates before calling cascade fit
400 std::vector<const xAOD::Vertex*> selectedPsi2Candidates;
401 for(auto vxcItr=psi2Container.cptr()->cbegin(); vxcItr!=psi2Container.cptr()->cend(); ++vxcItr) {
402 // Check the passed flag first
403 const xAOD::Vertex* vtx = *vxcItr;
404 bool passed = false;
405 for(size_t i=0; i<m_vertexPsi2HypoNames.size(); i++) {
406 SG::AuxElement::Accessor<Char_t> flagAcc("passed_"+m_vertexPsi2HypoNames[i]);
407 if(flagAcc.isAvailable(*vtx) && flagAcc(*vtx)) {
408 passed |= 1;
409 }
410 }
411 if(m_vertexPsi2HypoNames.size() && !passed) continue;
412
413 // Check Psi2 candidate invariant mass and skip if need be
414 double mass_psi2 = m_V0Tools->invariantMass(*vxcItr, massesPsi2);
415 if (mass_psi2 < m_psi2MassLower || mass_psi2 > m_psi2MassUpper) continue;
416
417 // Add loose cut on Jpsi2 mass from Psi2 -> Jpsi2 pi+ pi-, or on phi mass from Ds+ -> phi pi+
418 TLorentzVector p4_mu1, p4_mu2;
419 p4_mu1.SetPtEtaPhiM( (*vxcItr)->trackParticle(0)->pt(),
420 (*vxcItr)->trackParticle(0)->eta(),
421 (*vxcItr)->trackParticle(0)->phi(), m_vtx2Daug1MassHypo);
422 p4_mu2.SetPtEtaPhiM( (*vxcItr)->trackParticle(1)->pt(),
423 (*vxcItr)->trackParticle(1)->eta(),
424 (*vxcItr)->trackParticle(1)->phi(), m_vtx2Daug2MassHypo);
425 double mass_jpsi2 = (p4_mu1 + p4_mu2).M();
426 if (mass_jpsi2 < m_jpsi2MassLower || mass_jpsi2 > m_jpsi2MassUpper) continue;
427
429 TLorentzVector p4_trk1, p4_trk2;
430 p4_trk1.SetPtEtaPhiM( (*vxcItr)->trackParticle(2)->pt(),
431 (*vxcItr)->trackParticle(2)->eta(),
432 (*vxcItr)->trackParticle(2)->phi(), m_vtx2Daug3MassHypo);
433 p4_trk2.SetPtEtaPhiM( (*vxcItr)->trackParticle(3)->pt(),
434 (*vxcItr)->trackParticle(3)->eta(),
435 (*vxcItr)->trackParticle(3)->phi(), m_vtx2Daug4MassHypo);
436 double mass_diTrk2 = (p4_trk1 + p4_trk2).M();
437 if (mass_diTrk2 < m_diTrack2MassLower || mass_diTrk2 > m_diTrack2MassUpper) continue;
438 }
439
440 double chi2DOF = (*vxcItr)->chiSquared()/(*vxcItr)->numberDoF();
441 if(m_chi2cut_Psi2>0 && chi2DOF>m_chi2cut_Psi2) continue;
442
443 selectedPsi2Candidates.push_back(*vxcItr);
444 }
445 if(selectedPsi2Candidates.size()==0) return StatusCode::SUCCESS;
446
447 // Select the Psi1 candidates before calling cascade fit
448 std::vector<const xAOD::Vertex*> selectedPsi1Candidates;
449 for(auto vxcItr=psi1Container.cptr()->cbegin(); vxcItr!=psi1Container.cptr()->cend(); ++vxcItr) {
450 // Check the passed flag first
451 const xAOD::Vertex* vtx = *vxcItr;
452 bool passed = false;
453 for(size_t i=0; i<m_vertexPsi1HypoNames.size(); i++) {
454 SG::AuxElement::Accessor<Char_t> flagAcc("passed_"+m_vertexPsi1HypoNames[i]);
455 if(flagAcc.isAvailable(*vtx) && flagAcc(*vtx)) {
456 passed |= 1;
457 }
458 }
459 if(m_vertexPsi1HypoNames.size() && !passed) continue;
460
461 // Check Psi candidate invariant mass and skip if need be
462 double mass_psi1 = m_V0Tools->invariantMass(*vxcItr,massesPsi1);
463 if(mass_psi1 < m_psi1MassLower || mass_psi1 > m_psi1MassUpper) continue;
464
465 // Add loose cut on Jpsi1 mass from Psi1 -> Jpsi1 pi+ pi-, or on phi mass from Ds+ -> phi pi+
466 TLorentzVector p4_mu1, p4_mu2;
467 p4_mu1.SetPtEtaPhiM( (*vxcItr)->trackParticle(0)->pt(),
468 (*vxcItr)->trackParticle(0)->eta(),
469 (*vxcItr)->trackParticle(0)->phi(), m_vtx1Daug1MassHypo);
470 p4_mu2.SetPtEtaPhiM( (*vxcItr)->trackParticle(1)->pt(),
471 (*vxcItr)->trackParticle(1)->eta(),
472 (*vxcItr)->trackParticle(1)->phi(), m_vtx1Daug2MassHypo);
473 double mass_jpsi1 = (p4_mu1 + p4_mu2).M();
474 if (mass_jpsi1 < m_jpsi1MassLower || mass_jpsi1 > m_jpsi1MassUpper) continue;
475
477 TLorentzVector p4_trk1, p4_trk2;
478 p4_trk1.SetPtEtaPhiM( (*vxcItr)->trackParticle(2)->pt(),
479 (*vxcItr)->trackParticle(2)->eta(),
480 (*vxcItr)->trackParticle(2)->phi(), m_vtx1Daug3MassHypo);
481 p4_trk2.SetPtEtaPhiM( (*vxcItr)->trackParticle(3)->pt(),
482 (*vxcItr)->trackParticle(3)->eta(),
483 (*vxcItr)->trackParticle(3)->phi(), m_vtx1Daug4MassHypo);
484 double mass_diTrk1 = (p4_trk1 + p4_trk2).M();
485 if (mass_diTrk1 < m_diTrack1MassLower || mass_diTrk1 > m_diTrack1MassUpper) continue;
486 }
487
488 double chi2DOF = (*vxcItr)->chiSquared()/(*vxcItr)->numberDoF();
489 if(m_chi2cut_Psi1>0 && chi2DOF>m_chi2cut_Psi1) continue;
490
491 selectedPsi1Candidates.push_back(*vxcItr);
492 }
493 if(selectedPsi1Candidates.size()==0) return StatusCode::SUCCESS;
494
495 std::vector<std::pair<const xAOD::Vertex*, const xAOD::Vertex*> > candidatePairs;
496 for(auto psi1Itr=selectedPsi1Candidates.cbegin(); psi1Itr!=selectedPsi1Candidates.cend(); ++psi1Itr) {
497 tracksPsi1.clear();
498 for(size_t i=0; i<(*psi1Itr)->nTrackParticles(); i++) tracksPsi1.push_back((*psi1Itr)->trackParticle(i));
499 for(auto psi2Itr=selectedPsi2Candidates.cbegin(); psi2Itr!=selectedPsi2Candidates.cend(); ++psi2Itr) {
500 bool skip = false;
501 for(size_t j=0; j<(*psi2Itr)->nTrackParticles(); j++) {
502 if(std::find(tracksPsi1.cbegin(), tracksPsi1.cend(), (*psi2Itr)->trackParticle(j)) != tracksPsi1.cend()) { skip = true; break; }
503 }
504 if(skip) continue;
506 for(size_t ic=0; ic<candidatePairs.size(); ic++) {
507 const xAOD::Vertex* psi1Vertex = candidatePairs[ic].first;
508 const xAOD::Vertex* psi2Vertex = candidatePairs[ic].second;
509 if((psi1Vertex == *psi1Itr && psi2Vertex == *psi2Itr) || (psi1Vertex == *psi2Itr && psi2Vertex == *psi1Itr)) { skip = true; break; }
510 }
511 }
512 if(skip) continue;
513 candidatePairs.push_back(std::pair<const xAOD::Vertex*, const xAOD::Vertex*>(*psi1Itr,*psi2Itr));
514 }
515 }
516
517 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(); } );
518 if(m_maxCandidates>0 && candidatePairs.size()>m_maxCandidates) {
519 candidatePairs.erase(candidatePairs.begin()+m_maxCandidates, candidatePairs.end());
520 }
521
522 for(size_t ic=0; ic<candidatePairs.size(); ic++) {
523 const xAOD::Vertex* psi1Vertex = candidatePairs[ic].first;
524 const xAOD::Vertex* psi2Vertex = candidatePairs[ic].second;
525
526 tracksPsi1.clear();
527 for(size_t it=0; it<psi1Vertex->nTrackParticles(); it++) tracksPsi1.push_back(psi1Vertex->trackParticle(it));
528 if (tracksPsi1.size() != massesPsi1.size()) {
529 ATH_MSG_ERROR("Problems with Psi1 input: number of tracks or track mass inputs is not correct!");
530 }
531 tracksPsi2.clear();
532 for(size_t it=0; it<psi2Vertex->nTrackParticles(); it++) tracksPsi2.push_back(psi2Vertex->trackParticle(it));
533 if (tracksPsi2.size() != massesPsi2.size()) {
534 ATH_MSG_ERROR("Problems with Psi2 input: number of tracks or track mass inputs is not correct!");
535 }
536
537 tracksJpsi1.clear();
538 tracksJpsi1.push_back(psi1Vertex->trackParticle(0));
539 tracksJpsi1.push_back(psi1Vertex->trackParticle(1));
540 tracksDiTrk1.clear();
541 if(m_vtx1Daug_num==4) {
542 tracksDiTrk1.push_back(psi1Vertex->trackParticle(2));
543 tracksDiTrk1.push_back(psi1Vertex->trackParticle(3));
544 }
545 tracksJpsi2.clear();
546 tracksJpsi2.push_back(psi2Vertex->trackParticle(0));
547 tracksJpsi2.push_back(psi2Vertex->trackParticle(1));
548 tracksDiTrk2.clear();
549 if(m_vtx2Daug_num==4) {
550 tracksDiTrk2.push_back(psi2Vertex->trackParticle(2));
551 tracksDiTrk2.push_back(psi2Vertex->trackParticle(3));
552 }
553
554 TLorentzVector p4_moth;
555 TLorentzVector tmp;
556 for(size_t it=0; it<psi1Vertex->nTrackParticles(); it++) {
557 tmp.SetPtEtaPhiM(psi1Vertex->trackParticle(it)->pt(),psi1Vertex->trackParticle(it)->eta(),psi1Vertex->trackParticle(it)->phi(),massesPsi1[it]);
558 p4_moth += tmp;
559 }
560 for(size_t it=0; it<psi2Vertex->nTrackParticles(); it++) {
561 tmp.SetPtEtaPhiM(psi2Vertex->trackParticle(it)->pt(),psi2Vertex->trackParticle(it)->eta(),psi2Vertex->trackParticle(it)->phi(),massesPsi2[it]);
562 p4_moth += tmp;
563 }
564 if (p4_moth.M() < m_MassLower || p4_moth.M() > m_MassUpper) continue;
565
566 // Apply the user's settings to the fitter
567 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
568 // Robustness: http://cdsweb.cern.ch/record/685551
569 int robustness = 0;
570 m_iVertexFitter->setRobustness(robustness, *state);
571 // Build up the topology
572 // Vertex list
573 std::vector<Trk::VertexID> vrtList;
574 // Psi1 vertex
575 Trk::VertexID vID1;
576 // https://gitlab.cern.ch/atlas/athena/-/blob/21.2/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
577 if (m_constrPsi1) {
578 vID1 = m_iVertexFitter->startVertex(tracksPsi1,massesPsi1,*state,m_massPsi1);
579 } else {
580 vID1 = m_iVertexFitter->startVertex(tracksPsi1,massesPsi1,*state);
581 }
582 vrtList.push_back(vID1);
583 // Psi2 vertex
584 Trk::VertexID vID2;
585 if (m_constrPsi2) {
586 vID2 = m_iVertexFitter->nextVertex(tracksPsi2,massesPsi2,*state,m_massPsi2);
587 } else {
588 vID2 = m_iVertexFitter->nextVertex(tracksPsi2,massesPsi2,*state);
589 }
590 vrtList.push_back(vID2);
591 // Mother vertex including Psi1 and Psi2
592 std::vector<const xAOD::TrackParticle*> tp; tp.clear();
593 std::vector<double> tp_masses; tp_masses.clear();
594 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList,*state);
595 if (m_constrJpsi1) {
596 std::vector<Trk::VertexID> cnstV; cnstV.clear();
597 if ( !m_iVertexFitter->addMassConstraint(vID1,tracksJpsi1,cnstV,*state,m_massJpsi1).isSuccess() ) {
598 ATH_MSG_WARNING("addMassConstraint for Jpsi1 failed");
599 }
600 }
602 std::vector<Trk::VertexID> cnstV; cnstV.clear();
603 if ( !m_iVertexFitter->addMassConstraint(vID1,tracksDiTrk1,cnstV,*state,m_massDiTrk1).isSuccess() ) {
604 ATH_MSG_WARNING("addMassConstraint for DiTrk1 failed");
605 }
606 }
607 if (m_constrJpsi2) {
608 std::vector<Trk::VertexID> cnstV; cnstV.clear();
609 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJpsi2,cnstV,*state,m_massJpsi2).isSuccess() ) {
610 ATH_MSG_WARNING("addMassConstraint for Jpsi2 failed");
611 }
612 }
614 std::vector<Trk::VertexID> cnstV; cnstV.clear();
615 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksDiTrk2,cnstV,*state,m_massDiTrk2).isSuccess() ) {
616 ATH_MSG_WARNING("addMassConstraint for DiTrk2 failed");
617 }
618 }
619 // Do the work
620 std::unique_ptr<Trk::VxCascadeInfo> result(m_iVertexFitter->fitCascade(*state));
621
622 bool pass = false;
623 if (result != nullptr) {
624 for(auto v : result->vertices()) {
625 if(v->nTrackParticles()==0) {
626 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
627 v->setTrackParticleLinks(nullLinkVector);
628 }
629 }
630 // reset links to original tracks
631 BPhysPVCascadeTools::PrepareVertexLinks(result.get(), trackContainer.cptr());
632
633 // necessary to prevent memory leak
634 result->setSVOwnership(true);
635
636 // Chi2/DOF cut
637 double chi2DOF = result->fitChi2()/result->nDoF();
638 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
639
640 if(chi2CutPassed) {
641 cascadeinfoContainer->push_back(result.release());
642 pass = true;
643 }
644 }
645
646 // do cascade fit again without any mass constraints
647 if(pass) {
649 std::unique_ptr<Trk::IVKalState> state (m_iVertexFitter->makeState(ctx));
650 m_iVertexFitter->setRobustness(robustness, *state);
651 std::vector<Trk::VertexID> vrtList_nc;
652 // Psi1 vertex
653 Trk::VertexID vID1_nc = m_iVertexFitter->startVertex(tracksPsi1,massesPsi1,*state);
654 vrtList_nc.push_back(vID1_nc);
655 // Psi2 vertex
656 Trk::VertexID vID2_nc = m_iVertexFitter->nextVertex(tracksPsi2,massesPsi2,*state);
657 vrtList_nc.push_back(vID2_nc);
658 // Mother vertex including Psi1 and Psi2
659 std::vector<const xAOD::TrackParticle*> tp; tp.clear();
660 std::vector<double> tp_masses; tp_masses.clear();
661 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList_nc,*state);
662 // Do the work
663 std::unique_ptr<Trk::VxCascadeInfo> result_nc(m_iVertexFitter->fitCascade(*state));
664
665 if (result_nc != nullptr) {
666 for(auto v : result_nc->vertices()) {
667 if(v->nTrackParticles()==0) {
668 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
669 v->setTrackParticleLinks(nullLinkVector);
670 }
671 }
672 // reset links to original tracks
673 BPhysPVCascadeTools::PrepareVertexLinks(result_nc.get(), trackContainer.cptr());
674
675 // necessary to prevent memory leak
676 result_nc->setSVOwnership(true);
677 cascadeinfoContainer_noConstr->push_back(result_nc.release());
678 }
679 else cascadeinfoContainer_noConstr->push_back(0);
680 }
681 else cascadeinfoContainer_noConstr->push_back(0);
682 }
683 } //Iterate over candidatePairs
684
685 return StatusCode::SUCCESS;
686 }
687}
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_ERROR(x)
#define ATH_MSG_FATAL(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.
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 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 const xAOD::Vertex * FindVertex(const xAOD::VertexContainer *c, const xAOD::Vertex *v)
static void SetVectorInfo(xAOD::BPhysHelper &, const Trk::VxCascadeInfo *)
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
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
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.