ATLAS Offline Software
Loading...
Searching...
No Matches
JpsiXPlusDisplaced.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 "JpsiXPlusDisplaced.h"
6#include "TrkVertexFitterInterfaces/IVertexFitter.h"
7#include "TrkVKalVrtFitter/TrkVKalVrtFitter.h"
8#include "TrkVKalVrtFitter/VxCascadeInfo.h"
9#include "TrkVertexAnalysisUtils/V0Tools.h"
10#include "TrkExInterfaces/IExtrapolator.h"
11#include "DerivationFrameworkBPhys/CascadeTools.h"
12#include "BPhysPVCascadeTools.h"
13#include "xAODTracking/VertexAuxContainer.h"
14#include "InDetConversionFinderTools/VertexPointEstimator.h"
15#include "xAODBPhys/BPhysHypoHelper.h"
16#include "HepPDT/ParticleDataTable.hh"
17#include "VxVertex/RecVertex.h"
18#include "JpsiUpsilonTools/JpsiUpsilonCommon.h"
19#include "TruthUtils/HepMCHelpers.h"
20#include <algorithm>
21#include <functional>
22
23namespace DerivationFramework {
24 typedef ElementLink<xAOD::VertexContainer> VertexLink;
25 typedef std::vector<VertexLink> VertexLinkVector;
26
27 using Analysis::JpsiUpsilonCommon;
28
29 JpsiXPlusDisplaced::MesonCandidateVector::MesonCandidateVector(size_t num, bool orderByPt = true) {
30 m_num = num;
31 m_orderByPt = orderByPt;
32 }
33
34 void JpsiXPlusDisplaced::MesonCandidateVector::push_back(const MesonCandidate& etac) {
35 if(m_num>0 && m_vector.size()>=m_num) {
36 if(m_orderByPt && etac.pt<=m_vector.back().pt) return;
37 else if(!m_orderByPt && etac.chi2NDF>=m_vector.back().chi2NDF) return;
38 }
39 auto pos = m_vector.cend();
40 for(auto iter=m_vector.cbegin(); iter!=m_vector.cend(); iter++) {
41 if(m_orderByPt) {
42 if(etac.pt>iter->pt) { pos = iter; break; }
43 }
44 else {
45 if(etac.chi2NDF<iter->chi2NDF) { pos = iter; break; }
46 }
47 }
48 m_vector.insert(pos, etac);
49 if(m_num>0 && m_vector.size()>m_num) m_vector.pop_back();
50 }
51
52 const std::vector<JpsiXPlusDisplaced::MesonCandidate>& JpsiXPlusDisplaced::MesonCandidateVector::vector() const {
53 return m_vector;
54 }
55
56 JpsiXPlusDisplaced::JpsiXPlusDisplaced(const std::string& type, const std::string& name, const IInterface* parent) : base_class(type,name,parent),
57 m_vertexJXContainerKey("InputJXVertices"),
58 m_vertexV0ContainerKey(""),
59 m_cascadeOutputKeys({"JpsiXPlusDisVtx1_sub", "JpsiXPlusDisVtx1", "JpsiXPlusDisVtx2", "JpsiXPlusDisVtx3"}),
60 m_cascadeOutputKeys_mvc({"JpsiXPlusDisVtx1_sub_mvc", "JpsiXPlusDisVtx1_mvc", "JpsiXPlusDisVtx2_mvc", "JpsiXPlusDisVtx3_mvc"}),
61 m_v0VtxOutputKey(""),
62 m_TrkParticleCollection("InDetTrackParticles"),
63 m_VxPrimaryCandidateName("PrimaryVertices"),
64 m_pvContainerName("PrimaryVertices"),
65 m_refPVContainerName("RefittedPrimaryVertices"),
66 m_eventInfo_key("EventInfo"),
67 m_RelinkContainers({"InDetTrackParticles","InDetLargeD0TrackParticles"}),
68 m_useImprovedMass(false),
69 m_jxMassLower(0.0),
70 m_jxMassUpper(10000.0),
71 m_jpsiMassLower(0.0),
72 m_jpsiMassUpper(10000.0),
73 m_diTrackMassLower(-1.0),
74 m_diTrackMassUpper(-1.0),
75 m_V0Hypothesis("Lambda"),
76 m_LambdaMassLower(0.0),
77 m_LambdaMassUpper(10000.0),
78 m_KsMassLower(0.0),
79 m_KsMassUpper(10000.0),
80 m_lxyV0_cut(-999.0),
81 m_minMass_gamma(-1.0),
82 m_chi2cut_gamma(-1.0),
83 m_DisplacedMassLower(0.0),
84 m_DisplacedMassUpper(10000.0),
85 m_lxyDisV_cut(-999.0),
86 m_lxyDpm_cut(-999.0),
87 m_lxyD0_cut(-999.0),
88 m_MassLower(0.0),
89 m_MassUpper(31000.0),
90 m_PostMassLower(0.0),
91 m_PostMassUpper(31000.0),
92 m_jxDaug_num(4),
93 m_jxDaug1MassHypo(-1),
94 m_jxDaug2MassHypo(-1),
95 m_jxDaug3MassHypo(-1),
96 m_jxDaug4MassHypo(-1),
97 m_jxPtOrdering(false),
98 m_disVDaug_num(3),
99 m_disVDaug3MassHypo(-1),
100 m_disVDaug3MinPt(480),
101 m_extraTrk1MassHypo(-1),
102 m_extraTrk1MinPt(480),
103 m_extraTrk2MassHypo(-1),
104 m_extraTrk2MinPt(480),
105 m_extraTrk3MassHypo(-1),
106 m_extraTrk3MinPt(480),
107 m_DpmMassLower(0.0),
108 m_DpmMassUpper(10000.0),
109 m_D0MassLower(0.0),
110 m_D0MassUpper(10000.0),
111 m_maxMesonCandidates(400),
112 m_MesonPtOrdering(true),
113 m_massJX(-1),
114 m_massJpsi(-1),
115 m_massX(-1),
116 m_massDisV(-1),
117 m_massLd(-1),
118 m_massKs(-1),
119 m_massDpm(-1),
120 m_massD0(-1),
121 m_massJXV0(-1),
122 m_massMainV(-1),
123 m_constrJX(false),
124 m_constrJpsi(false),
125 m_constrX(false),
126 m_constrDisV(false),
127 m_constrV0(false),
128 m_constrDpm(false),
129 m_constrD0(false),
130 m_constrJXV0(false),
131 m_constrMainV(false),
132 m_cascadeFitWithPV(0),
133 m_firstDecayAtPV(false),
134 m_doPostMainVContrFit(false),
135 m_JXSubVtx(false),
136 m_JXV0SubVtx(false),
137 m_chi2cut_JX(-1.0),
138 m_chi2cut_V0(-1.0),
139 m_chi2cut_DisV(-1.0),
140 m_chi2cut_Dpm(-1.0),
141 m_chi2cut_D0(-1.0),
142 m_chi2cut(-1.0),
143 m_useTRT(false),
144 m_ptTRT(450),
145 m_d0_cut(2),
146 m_maxJXCandidates(0),
147 m_maxV0Candidates(0),
148 m_maxDisVCandidates(0),
149 m_maxMainVCandidates(0),
150 m_iVertexFitter("Trk::TrkVKalVrtFitter"),
151 m_iV0Fitter("Trk::V0VertexFitter"),
152 m_iGammaFitter("Trk::TrkVKalVrtFitter"),
153 m_pvRefitter("Analysis::PrimaryVertexRefitter"),
154 m_V0Tools("Trk::V0Tools"),
155 m_trackToVertexTool("Reco::TrackToVertex"),
156 m_trkSelector("InDet::TrackSelectorTool"),
157 m_v0TrkSelector("InDet::TrackSelectorTool"),
158 m_CascadeTools("DerivationFramework::CascadeTools"),
159 m_vertexEstimator("InDet::VertexPointEstimator"),
160 m_extrapolator("Trk::Extrapolator/AtlasExtrapolator")
161 {
162 declareProperty("JXVertices", m_vertexJXContainerKey);
163 declareProperty("V0Vertices", m_vertexV0ContainerKey);
164 declareProperty("JXVtxHypoNames", m_vertexJXHypoNames);
165 declareProperty("CascadeVertexCollections", m_cascadeOutputKeys); // size is 3 or 4 only
166 declareProperty("CascadeVertexCollectionsMVC",m_cascadeOutputKeys_mvc); // size is 3 or 4 only
167 declareProperty("OutputV0VtxCollection", m_v0VtxOutputKey);
168 declareProperty("TrackParticleCollection", m_TrkParticleCollection);
169 declareProperty("VxPrimaryCandidateName", m_VxPrimaryCandidateName);
170 declareProperty("PVContainerName", m_pvContainerName);
171 declareProperty("RefPVContainerName", m_refPVContainerName);
172 declareProperty("EventInfoKey", m_eventInfo_key);
173 declareProperty("RelinkTracks", m_RelinkContainers);
174 declareProperty("UseImprovedMass", m_useImprovedMass);
175 declareProperty("JXMassLowerCut", m_jxMassLower); // only effective when m_jxDaug_num>2
176 declareProperty("JXMassUpperCut", m_jxMassUpper); // only effective when m_jxDaug_num>2
177 declareProperty("JpsiMassLowerCut", m_jpsiMassLower);
178 declareProperty("JpsiMassUpperCut", m_jpsiMassUpper);
179 declareProperty("DiTrackMassLower", m_diTrackMassLower); // only effective when m_jxDaug_num=4
180 declareProperty("DiTrackMassUpper", m_diTrackMassUpper); // only effective when m_jxDaug_num=4
181 declareProperty("V0Hypothesis", m_V0Hypothesis); // "Ks" or "Lambda"
182 declareProperty("LambdaMassLowerCut", m_LambdaMassLower);
183 declareProperty("LambdaMassUpperCut", m_LambdaMassUpper);
184 declareProperty("KsMassLowerCut", m_KsMassLower);
185 declareProperty("KsMassUpperCut", m_KsMassUpper);
186 declareProperty("LxyV0Cut", m_lxyV0_cut);
187 declareProperty("MassCutGamma", m_minMass_gamma);
188 declareProperty("Chi2CutGamma", m_chi2cut_gamma);
189 declareProperty("DisplacedMassLowerCut", m_DisplacedMassLower); // only effective when m_disVDaug_num=3
190 declareProperty("DisplacedMassUpperCut", m_DisplacedMassUpper); // only effective when m_disVDaug_num=3
191 declareProperty("LxyDisVtxCut", m_lxyDisV_cut); // only effective when m_disVDaug_num=3
192 declareProperty("LxyDpmCut", m_lxyDpm_cut); // only effective for D+/-
193 declareProperty("LxyD0Cut", m_lxyD0_cut); // only effective for D0
194 declareProperty("MassLowerCut", m_MassLower);
195 declareProperty("MassUpperCut", m_MassUpper);
196 declareProperty("PostMassLowerCut", m_PostMassLower); // only effective when m_doPostMainVContrFit=true
197 declareProperty("PostMassUpperCut", m_PostMassUpper); // only effective when m_doPostMainVContrFit=true
198 declareProperty("HypothesisName", m_hypoName = "TQ");
199 declareProperty("NumberOfJXDaughters", m_jxDaug_num); // 2, or 3, or 4 only
200 declareProperty("JXDaug1MassHypo", m_jxDaug1MassHypo);
201 declareProperty("JXDaug2MassHypo", m_jxDaug2MassHypo);
202 declareProperty("JXDaug3MassHypo", m_jxDaug3MassHypo);
203 declareProperty("JXDaug4MassHypo", m_jxDaug4MassHypo);
204 declareProperty("JXPtOrdering", m_jxPtOrdering);
205 declareProperty("NumberOfDisVDaughters", m_disVDaug_num); // 2 or 3 only
206 declareProperty("DisVDaug3MassHypo", m_disVDaug3MassHypo); // only effective when m_disVDaug_num=3
207 declareProperty("DisVDaug3MinPt", m_disVDaug3MinPt); // only effective when m_disVDaug_num=3
208 declareProperty("ExtraTrack1MassHypo", m_extraTrk1MassHypo); // for decays like B- -> J/psi Lambda pbar, the extra track is pbar (for m_disVDaug_num=2 only now)
209 declareProperty("ExtraTrack1MinPt", m_extraTrk1MinPt); // only effective if m_extraTrk1MassHypo>0
210 declareProperty("ExtraTrack2MassHypo", m_extraTrk2MassHypo); // for decays like Xi_bc^0 -> Jpsi Lambda D0(->Kpi) (for m_disVDaug_num=2 only now)
211 declareProperty("ExtraTrack2MinPt", m_extraTrk2MinPt); // only effective if m_extraTrk2MassHypo>0
212 declareProperty("ExtraTrack3MassHypo", m_extraTrk3MassHypo); // for decays like Bc+ -> Jpsi D+ Ks (for m_disVDaug_num=2 only now)
213 declareProperty("ExtraTrack3MinPt", m_extraTrk3MinPt); // only effective if m_extraTrk3MassHypo>0
214 declareProperty("DpmMassLowerCut", m_DpmMassLower); // only for D+/-
215 declareProperty("DpmMassUpperCut", m_DpmMassUpper); // only for D+/-
216 declareProperty("D0MassLowerCut", m_D0MassLower); // only for D0
217 declareProperty("D0MassUpperCut", m_D0MassUpper); // only for D0
218 declareProperty("MaxMesonCandidates", m_maxMesonCandidates); // only for 2/3 extra tracks
219 declareProperty("MesonPtOrdering", m_MesonPtOrdering); // only for 2/3 extra tracks
220 declareProperty("JXMass", m_massJX); // only effective when m_jxDaug_num>2
221 declareProperty("JpsiMass", m_massJpsi);
222 declareProperty("XMass", m_massX); // only effective when m_jxDaug_num=4
223 declareProperty("DisVtxMass", m_massDisV); // only effective when m_disVDaug_num=3
224 declareProperty("LambdaMass", m_massLd);
225 declareProperty("KsMass", m_massKs);
226 declareProperty("DpmMass", m_massDpm);
227 declareProperty("D0Mass", m_massD0);
228 declareProperty("JXV0Mass", m_massJXV0);
229 declareProperty("MainVtxMass", m_massMainV);
230 declareProperty("ApplyJXMassConstraint", m_constrJX); // only effective when m_jxDaug_num>2
231 declareProperty("ApplyJpsiMassConstraint", m_constrJpsi);
232 declareProperty("ApplyXMassConstraint", m_constrX); // only effective when m_jxDaug_num=4
233 declareProperty("ApplyDisVMassConstraint", m_constrDisV); // only effective when m_disVDaug_num=3
234 declareProperty("ApplyV0MassConstraint", m_constrV0);
235 declareProperty("ApplyDpmMassConstraint", m_constrDpm); // only for D+/-
236 declareProperty("ApplyD0MassConstraint", m_constrD0); // only for D0
237 declareProperty("ApplyJXV0MassConstraint", m_constrJXV0);
238 declareProperty("ApplyMainVMassConstraint", m_constrMainV);
239 declareProperty("DoCascadeFitWithPV", m_cascadeFitWithPV);
240 declareProperty("FirstDecayAtPV", m_firstDecayAtPV);
241 declareProperty("DoPostMainVContrFit", m_doPostMainVContrFit); // only effective when m_constrMainV=false
242 declareProperty("HasJXSubVertex", m_JXSubVtx);
243 declareProperty("HasJXV0SubVertex", m_JXV0SubVtx);
244 declareProperty("Chi2CutJX", m_chi2cut_JX);
245 declareProperty("Chi2CutV0", m_chi2cut_V0);
246 declareProperty("Chi2CutDisV", m_chi2cut_DisV); // only effective when m_disVDaug_num=3
247 declareProperty("Chi2CutDpm", m_chi2cut_Dpm); // only for D+/-
248 declareProperty("Chi2CutD0", m_chi2cut_D0); // only for D0
249 declareProperty("Chi2Cut", m_chi2cut);
250 declareProperty("UseTRT", m_useTRT);
251 declareProperty("PtTRT", m_ptTRT);
252 declareProperty("Trackd0Cut", m_d0_cut);
253 declareProperty("MaxJXCandidates", m_maxJXCandidates);
254 declareProperty("MaxV0Candidates", m_maxV0Candidates);
255 declareProperty("MaxDisVCandidates", m_maxDisVCandidates); // only effective when m_disVDaug_num=3
256 declareProperty("MaxMainVCandidates", m_maxMainVCandidates);
257 declareProperty("RefitPV", m_refitPV = true);
258 declareProperty("MaxnPV", m_PV_max = 1000);
259 declareProperty("MinNTracksInPV", m_PV_minNTracks = 0);
260 declareProperty("DoVertexType", m_DoVertexType = 7);
261 declareProperty("TrkVertexFitterTool", m_iVertexFitter);
262 declareProperty("V0VertexFitterTool", m_iV0Fitter);
263 declareProperty("GammaFitterTool", m_iGammaFitter);
264 declareProperty("PVRefitter", m_pvRefitter);
265 declareProperty("V0Tools", m_V0Tools);
266 declareProperty("TrackToVertexTool", m_trackToVertexTool);
267 declareProperty("TrackSelectorTool", m_trkSelector);
268 declareProperty("V0TrackSelectorTool", m_v0TrkSelector);
269 declareProperty("CascadeTools", m_CascadeTools);
270 declareProperty("VertexPointEstimator", m_vertexEstimator);
271 declareProperty("Extrapolator", m_extrapolator);
272 }
273
274 StatusCode JpsiXPlusDisplaced::initialize() {
275 if(m_V0Hypothesis != "Ks" && m_V0Hypothesis != "Lambda"
276 && m_V0Hypothesis != "Lambda/Ks"&& m_V0Hypothesis != "Ks/Lambda") {
277 ATH_MSG_FATAL("Incorrect V0 container hypothesis - not recognized");
278 return StatusCode::FAILURE;
279 }
280
281 if(m_jxDaug_num<2 || m_jxDaug_num>4 || m_disVDaug_num<2 || m_disVDaug_num>3) {
282 ATH_MSG_FATAL("Incorrect number of JX or DisVtx daughters");
283 return StatusCode::FAILURE;
284 }
285
286 if(m_vertexV0ContainerKey.key()=="" && m_v0VtxOutputKey.key()=="") {
287 ATH_MSG_FATAL("Input and output V0 container names can not be both empty");
288 return StatusCode::FAILURE;
289 }
290
291 // retrieving vertex Fitter
292 ATH_CHECK( m_iVertexFitter.retrieve() );
293
294 // retrieving V0 vertex Fitter
295 ATH_CHECK( m_iV0Fitter.retrieve() );
296
297 // retrieving photon conversion vertex Fitter
298 ATH_CHECK( m_iGammaFitter.retrieve() );
299
300 // retrieving primary vertex refitter
301 ATH_CHECK( m_pvRefitter.retrieve() );
302
303 // retrieving the V0 tool
304 ATH_CHECK( m_V0Tools.retrieve() );
305
306 // retrieving the TrackToVertex extrapolator tool
307 ATH_CHECK( m_trackToVertexTool.retrieve() );
308
309 // retrieving the track selector tool
310 ATH_CHECK( m_trkSelector.retrieve() );
311
312 // retrieving the V0 track selector tool
313 ATH_CHECK( m_v0TrkSelector.retrieve() );
314
315 // retrieving the Cascade tools
316 ATH_CHECK( m_CascadeTools.retrieve() );
317
318 // retrieving the vertex point estimator
319 ATH_CHECK( m_vertexEstimator.retrieve() );
320
321 // retrieving the extrapolator
322 ATH_CHECK( m_extrapolator.retrieve() );
323
324 ATH_CHECK( m_vertexJXContainerKey.initialize() );
325 ATH_CHECK( m_vertexV0ContainerKey.initialize(SG::AllowEmpty) );
326 ATH_CHECK( m_VxPrimaryCandidateName.initialize() );
327 ATH_CHECK( m_pvContainerName.initialize() );
328 ATH_CHECK( m_refPVContainerName.initialize() );
329 ATH_CHECK( m_TrkParticleCollection.initialize() );
330 ATH_CHECK( m_cascadeOutputKeys.initialize() );
331 ATH_CHECK( m_cascadeOutputKeys_mvc.initialize() );
332 ATH_CHECK( m_eventInfo_key.initialize() );
333 ATH_CHECK( m_RelinkContainers.initialize() );
334 ATH_CHECK( m_v0VtxOutputKey.initialize(SG::AllowEmpty) );
335
336 ATH_CHECK( m_partPropSvc.retrieve() );
337 auto pdt = m_partPropSvc->PDT();
338
339 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Generators/TruthUtils/TruthUtils/AtlasPID.h
340 m_mass_e = BPhysPVCascadeTools::getParticleMass(pdt, MC::ELECTRON);
341 m_mass_mu = BPhysPVCascadeTools::getParticleMass(pdt, MC::MUON);
342 m_mass_pion = BPhysPVCascadeTools::getParticleMass(pdt, MC::PIPLUS);
343 m_mass_proton = BPhysPVCascadeTools::getParticleMass(pdt, MC::PROTON);
344 m_mass_Lambda = BPhysPVCascadeTools::getParticleMass(pdt, MC::LAMBDA0);
345 m_mass_Ks = BPhysPVCascadeTools::getParticleMass(pdt, MC::K0S);
346 m_mass_Xi = BPhysPVCascadeTools::getParticleMass(pdt, 3312);
347 m_mass_phi = BPhysPVCascadeTools::getParticleMass(pdt, 333);
348 m_mass_B0 = BPhysPVCascadeTools::getParticleMass(pdt, MC::B0);
349 m_mass_Dpm = BPhysPVCascadeTools::getParticleMass(pdt, MC::DPLUS);
350 m_mass_D0 = BPhysPVCascadeTools::getParticleMass(pdt, MC::D0);
351 m_mass_BCPLUS = BPhysPVCascadeTools::getParticleMass(pdt, MC::BCPLUS);
352 m_mass_Lambdab = BPhysPVCascadeTools::getParticleMass(pdt, MC::LAMBDAB0);
353
354 m_massesV0_ppi.push_back(m_mass_proton);
355 m_massesV0_ppi.push_back(m_mass_pion);
356 m_massesV0_pip.push_back(m_mass_pion);
357 m_massesV0_pip.push_back(m_mass_proton);
358 m_massesV0_pipi.push_back(m_mass_pion);
359 m_massesV0_pipi.push_back(m_mass_pion);
360
361 // retrieve particle masses
362 if(m_constrJpsi && m_massJpsi<0) m_massJpsi = BPhysPVCascadeTools::getParticleMass(pdt, MC::JPSI);
363 if(m_jxDaug_num>=3 && m_constrJX && m_massJX<0) m_massJX = BPhysPVCascadeTools::getParticleMass(pdt, MC::PSI2S);
364 if(m_jxDaug_num==4 && m_constrX && m_massX<0) m_massX = m_mass_phi;
365 if(m_constrV0) {
366 if(m_massLd<0) m_massLd = m_mass_Lambda;
367 if(m_massKs<0) m_massKs = m_mass_Ks;
368 }
369 if(m_disVDaug_num==3 && m_constrDisV && m_massDisV<0) m_massDisV = m_mass_Xi;
370 if(m_massMainV<0) m_massMainV = m_mass_BCPLUS;
371 if(m_constrDpm && m_massDpm<0) m_massDpm = m_mass_Dpm;
372 if(m_constrD0 && m_massD0<0) m_massD0 = m_mass_D0;
373 if(m_constrJXV0 && m_massJXV0<0) m_massJXV0 = m_mass_Lambdab;
374
375 if(m_jxDaug1MassHypo < 0.) m_jxDaug1MassHypo = m_mass_mu;
376 if(m_jxDaug2MassHypo < 0.) m_jxDaug2MassHypo = m_mass_mu;
377 if(m_jxDaug_num>=3 && m_jxDaug3MassHypo < 0.) m_jxDaug3MassHypo = m_mass_pion;
378 if(m_jxDaug_num==4 && m_jxDaug4MassHypo < 0.) m_jxDaug4MassHypo = m_mass_pion;
379 if(m_disVDaug_num==3 && m_disVDaug3MassHypo < 0.) m_disVDaug3MassHypo = m_mass_pion;
380
381 return StatusCode::SUCCESS;
382 }
383
384 StatusCode JpsiXPlusDisplaced::performSearch(std::vector<std::pair<Trk::VxCascadeInfo*,Trk::VxCascadeInfo*> >& cascadeinfoContainer, const std::vector<std::pair<const xAOD::Vertex*,V0Enum> >& selectedV0Candidates, const std::vector<const xAOD::TrackParticle*>& tracksDisplaced, const EventContext& ctx) const {
385 ATH_MSG_DEBUG( "JpsiXPlusDisplaced::performSearch" );
386 if(selectedV0Candidates.size()==0) return StatusCode::SUCCESS;
387
388 // Get TrackParticle container (standard + LRT)
389 SG::ReadHandle<xAOD::TrackParticleContainer> trackContainer(m_TrkParticleCollection, ctx);
390 ATH_CHECK( trackContainer.isValid() );
391
392 // Get all track containers when m_RelinkContainers is not empty
393 std::vector<const xAOD::TrackParticleContainer*> trackCols;
394 for(const SG::ReadHandleKey<xAOD::TrackParticleContainer>& key : m_RelinkContainers){
395 SG::ReadHandle<xAOD::TrackParticleContainer> handle(key, ctx);
396 ATH_CHECK( handle.isValid() );
397 trackCols.push_back(handle.cptr());
398 }
399
400 // Get default PV container
401 SG::ReadHandle<xAOD::VertexContainer> defaultPVContainer(m_VxPrimaryCandidateName, ctx);
402 ATH_CHECK( defaultPVContainer.isValid() );
403
404 // Get PV container
405 SG::ReadHandle<xAOD::VertexContainer> pvContainer(m_pvContainerName,ctx);
406 ATH_CHECK( pvContainer.isValid() );
407
408 // Get Jpsi+X container
409 SG::ReadHandle<xAOD::VertexContainer> jxContainer(m_vertexJXContainerKey, ctx);
410 ATH_CHECK( jxContainer.isValid() );
411
412 std::vector<double> massesJX{m_jxDaug1MassHypo, m_jxDaug2MassHypo};
413 if(m_jxDaug_num>=3) massesJX.push_back(m_jxDaug3MassHypo);
414 if(m_jxDaug_num==4) massesJX.push_back(m_jxDaug4MassHypo);
415
416 // Make the displaced candidates if needed
417 std::vector<XiCandidate> disVtxContainer;
418 if(m_disVDaug_num==3) {
419 for(size_t it=0; it<selectedV0Candidates.size(); ++it) {
420 std::pair<const xAOD::Vertex*,V0Enum> elem = selectedV0Candidates[it];
421 std::vector<const xAOD::TrackParticle*> tracksV0;
422 tracksV0.reserve(elem.first->nTrackParticles());
423 for(size_t i=0; i<elem.first->nTrackParticles(); i++) tracksV0.push_back(elem.first->trackParticle(i));
424 std::vector<double> massesV0;
425 if(elem.second==LAMBDA) massesV0 = m_massesV0_ppi;
426 else if(elem.second==LAMBDABAR) massesV0 = m_massesV0_pip;
427 else if(elem.second==KS) massesV0 = m_massesV0_pipi;
428 xAOD::BPhysHelper V0_helper(elem.first); TLorentzVector p4_v0;
429 for(int i=0; i<V0_helper.nRefTrks(); i++) p4_v0 += V0_helper.refTrk(i,massesV0[i]);
430 for(const xAOD::TrackParticle* TP : tracksDisplaced) {
431 if(TP->pt() < m_disVDaug3MinPt) continue;
432 // Check overlap
433 if(std::find(tracksV0.cbegin(), tracksV0.cend(), TP) != tracksV0.cend()) continue;
434 TLorentzVector tmp;
435 tmp.SetPtEtaPhiM(TP->pt(),TP->eta(),TP->phi(),m_disVDaug3MassHypo);
436 double disV_mass = (p4_v0+tmp).M();
437 if(m_useImprovedMass) {
438 if((elem.second==LAMBDA || elem.second==LAMBDABAR) && m_massLd>0) disV_mass += - p4_v0.M() + m_massLd;
439 else if(elem.second==KS && m_massKs>0) disV_mass += - p4_v0.M() + m_massKs;
440 }
441 // A rough mass window cut, as V0 and track3 are not from a common vertex
442 if(disV_mass > m_DisplacedMassLower-500. && disV_mass < m_DisplacedMassUpper+500.) {
443 auto disVtx = getXiCandidate(ctx,elem.first,elem.second,TP);
444 if(disVtx.V0vtx && disVtx.track) disVtxContainer.push_back(disVtx);
445 }
446 }
447 }
448
449 std::sort( disVtxContainer.begin(), disVtxContainer.end(), [](const XiCandidate& a, const XiCandidate& b) { return a.chi2NDF < b.chi2NDF; } );
450 if(m_maxDisVCandidates>0 && disVtxContainer.size()>m_maxDisVCandidates) {
451 disVtxContainer.erase(disVtxContainer.begin()+m_maxDisVCandidates, disVtxContainer.end());
452 }
453 if(disVtxContainer.size()==0) return StatusCode::SUCCESS;
454 } // m_disVDaug_num==3
455
456 // Select the JX candidates before calling cascade fit
457 std::vector<const xAOD::Vertex*> selectedJXCandidates;
458 for(const xAOD::Vertex* vtx : *jxContainer.cptr()) {
459 // Check the passed flag first
460 bool passed = false;
461 for(const std::string& name : m_vertexJXHypoNames) {
462 SG::AuxElement::Accessor<Char_t> flagAcc("passed_"+name);
463 if(flagAcc.isAvailable(*vtx) && flagAcc(*vtx)) {
464 passed = true;
465 }
466 }
467 if(m_vertexJXHypoNames.size() && !passed) continue;
468
469 // Add loose cut on Jpsi mass from e.g. JX -> Jpsi pi+ pi-
470 TLorentzVector p4_mu1, p4_mu2;
471 p4_mu1.SetPtEtaPhiM(vtx->trackParticle(0)->pt(),vtx->trackParticle(0)->eta(),vtx->trackParticle(0)->phi(), m_jxDaug1MassHypo);
472 p4_mu2.SetPtEtaPhiM(vtx->trackParticle(1)->pt(),vtx->trackParticle(1)->eta(),vtx->trackParticle(1)->phi(), m_jxDaug2MassHypo);
473 double mass_jpsi = (p4_mu1 + p4_mu2).M();
474 if (mass_jpsi < m_jpsiMassLower || mass_jpsi > m_jpsiMassUpper) continue;
475
476 TLorentzVector p4_trk1, p4_trk2;
477 if(m_jxDaug_num>=3) p4_trk1.SetPtEtaPhiM(vtx->trackParticle(2)->pt(),vtx->trackParticle(2)->eta(),vtx->trackParticle(2)->phi(), m_jxDaug3MassHypo);
478 if(m_jxDaug_num==4) p4_trk2.SetPtEtaPhiM(vtx->trackParticle(3)->pt(),vtx->trackParticle(3)->eta(),vtx->trackParticle(3)->phi(), m_jxDaug4MassHypo);
479
480 if(m_jxDaug_num==3) {
481 double mass_jx = (p4_mu1 + p4_mu2 + p4_trk1).M();
482 if(m_useImprovedMass && m_massJpsi>0) mass_jx += - (p4_mu1 + p4_mu2).M() + m_massJpsi;
483 if(mass_jx < m_jxMassLower || mass_jx > m_jxMassUpper) continue;
484 }
485 else if(m_jxDaug_num==4) {
486 double mass_jx = (p4_mu1 + p4_mu2 + p4_trk1 + p4_trk2).M();
487 if(m_useImprovedMass && m_massJpsi>0) mass_jx += - (p4_mu1 + p4_mu2).M() + m_massJpsi;
488 if(mass_jx < m_jxMassLower || mass_jx > m_jxMassUpper) continue;
489
490 if(m_diTrackMassLower>=0 && m_diTrackMassUpper>m_diTrackMassLower) {
491 double mass_diTrk = (p4_trk1 + p4_trk2).M();
492 if(mass_diTrk < m_diTrackMassLower || mass_diTrk > m_diTrackMassUpper) continue;
493 }
494 }
495
496 double chi2DOF = vtx->chiSquared()/vtx->numberDoF();
497 if(m_chi2cut_JX>0 && chi2DOF>m_chi2cut_JX) continue;
498
499 selectedJXCandidates.push_back(vtx);
500 }
501 if(selectedJXCandidates.size()==0) return StatusCode::SUCCESS;
502
503 if(m_jxPtOrdering) {
504 std::sort( selectedJXCandidates.begin(), selectedJXCandidates.end(), [massesJX](const xAOD::Vertex* a, const xAOD::Vertex* b) {
505 TLorentzVector p4_a, p4_b, tmp;
506 for(size_t it=0; it<a->nTrackParticles(); it++) {
507 tmp.SetPtEtaPhiM(a->trackParticle(it)->pt(), a->trackParticle(it)->eta(), a->trackParticle(it)->phi(), massesJX[it]);
508 p4_a += tmp;
509 }
510 for(size_t it=0; it<b->nTrackParticles(); it++) {
511 tmp.SetPtEtaPhiM(b->trackParticle(it)->pt(), b->trackParticle(it)->eta(), b->trackParticle(it)->phi(), massesJX[it]);
512 p4_b += tmp;
513 }
514 return p4_a.Pt() > p4_b.Pt();
515 } );
516 }
517 else {
518 std::sort( selectedJXCandidates.begin(), selectedJXCandidates.end(), [](const xAOD::Vertex* a, const xAOD::Vertex* b) { return a->chiSquared()/a->numberDoF() < b->chiSquared()/b->numberDoF(); } );
519 }
520 if(m_maxJXCandidates>0 && selectedJXCandidates.size()>m_maxJXCandidates) {
521 selectedJXCandidates.erase(selectedJXCandidates.begin()+m_maxJXCandidates, selectedJXCandidates.end());
522 }
523
524 // Select JX+DisV candidates
525 // Iterate over JX vertices
526 for(const xAOD::Vertex* jxVtx : selectedJXCandidates) {
527 // Iterate over displaced vertices
528 if(m_disVDaug_num==2) {
529 for(auto&& V0Candidate : selectedV0Candidates) {
530 std::vector<std::pair<Trk::VxCascadeInfo*, Trk::VxCascadeInfo*> > result = fitMainVtx(ctx, jxVtx, massesJX, V0Candidate.first, V0Candidate.second, trackContainer.cptr(), trackCols, defaultPVContainer.cptr(), pvContainer.cptr());
531 for(auto cascade_info_pair : result) {
532 if(cascade_info_pair.first) cascadeinfoContainer.push_back(cascade_info_pair);
533 }
534 }
535 } // m_disVDaug_num==2
536 else if(m_disVDaug_num==3) {
537 for(auto&& disVtx : disVtxContainer) {
538 std::vector<std::pair<Trk::VxCascadeInfo*, Trk::VxCascadeInfo*> > result = fitMainVtx(ctx, jxVtx, massesJX, disVtx, trackContainer.cptr(), trackCols, defaultPVContainer.cptr(), pvContainer.cptr());
539 for(auto cascade_info_pair : result) {
540 if(cascade_info_pair.first) cascadeinfoContainer.push_back(cascade_info_pair);
541 }
542 }
543 } // m_disVDaug_num==3
544 } // Iterate over JX vertices
545
546 return StatusCode::SUCCESS;
547 }
548
549 StatusCode JpsiXPlusDisplaced::addBranches(const EventContext& ctx) const {
550 size_t topoN = (m_disVDaug_num==2 ? 3 : 4);
551 if(!m_JXSubVtx) topoN--;
552 if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0) { // special cases
553 if(m_JXV0SubVtx) topoN = 4;
554 else topoN = 3;
555 }
556 else if(m_extraTrk1MassHypo>0 && m_disVDaug_num==2) { // special cases
557 if(m_JXV0SubVtx || m_JXSubVtx) topoN = 3;
558 else topoN = 2;
559 }
560
561 if(m_cascadeOutputKeys.size() != topoN) {
562 ATH_MSG_FATAL("Incorrect number of output cascade vertices");
563 return StatusCode::FAILURE;
564 }
565 if(!m_constrMainV && m_doPostMainVContrFit && m_cascadeOutputKeys_mvc.size() != topoN) {
566 ATH_MSG_FATAL("Incorrect number of output (mvc) cascade vertices");
567 return StatusCode::FAILURE;
568 }
569
570 std::array<SG::WriteHandle<xAOD::VertexContainer>, 4> VtxWriteHandles; int ikey(0);
571 for(const SG::WriteHandleKey<xAOD::VertexContainer>& key : m_cascadeOutputKeys) {
572 VtxWriteHandles[ikey] = SG::WriteHandle<xAOD::VertexContainer>(key, ctx);
573 ATH_CHECK( VtxWriteHandles[ikey].record(std::make_unique<xAOD::VertexContainer>(), std::make_unique<xAOD::VertexAuxContainer>()) );
574 ikey++;
575 }
576 std::array<SG::WriteHandle<xAOD::VertexContainer>, 4> VtxWriteHandles_mvc; int ikey_mvc(0);
577 if(!m_constrMainV && m_doPostMainVContrFit) {
578 for(const SG::WriteHandleKey<xAOD::VertexContainer>& key_mvc : m_cascadeOutputKeys_mvc) {
579 VtxWriteHandles_mvc[ikey_mvc] = SG::WriteHandle<xAOD::VertexContainer>(key_mvc, ctx);
580 ATH_CHECK( VtxWriteHandles_mvc[ikey_mvc].record(std::make_unique<xAOD::VertexContainer>(), std::make_unique<xAOD::VertexAuxContainer>()) );
581 ikey_mvc++;
582 }
583 }
584
585 //----------------------------------------------------
586 // retrieve primary vertices
587 //----------------------------------------------------
588 const xAOD::Vertex* primaryVertex(nullptr);
589 SG::ReadHandle<xAOD::VertexContainer> defaultPVContainer(m_VxPrimaryCandidateName, ctx);
590 ATH_CHECK( defaultPVContainer.isValid() );
591 if (defaultPVContainer.cptr()->size()==0) {
592 ATH_MSG_WARNING("You have no primary vertices: " << defaultPVContainer.cptr()->size());
593 return StatusCode::RECOVERABLE;
594 }
595 else primaryVertex = (*defaultPVContainer.cptr())[0];
596
597 //----------------------------------------------------
598 // Record refitted primary vertices
599 //----------------------------------------------------
600 SG::WriteHandle<xAOD::VertexContainer> refPvContainer;
601 if(m_refitPV) {
602 refPvContainer = SG::WriteHandle<xAOD::VertexContainer>(m_refPVContainerName, ctx);
603 ATH_CHECK( refPvContainer.record(std::make_unique<xAOD::VertexContainer>(), std::make_unique<xAOD::VertexAuxContainer>()) );
604 }
605
606 // Get TrackParticle container (standard + LRT)
607 SG::ReadHandle<xAOD::TrackParticleContainer> trackContainer(m_TrkParticleCollection, ctx);
608 ATH_CHECK( trackContainer.isValid() );
609
610 // Get all track containers when m_RelinkContainers is not empty
611 std::vector<const xAOD::TrackParticleContainer*> trackCols;
612 for(const SG::ReadHandleKey<xAOD::TrackParticleContainer>& key : m_RelinkContainers){
613 SG::ReadHandle<xAOD::TrackParticleContainer> handle(key, ctx);
614 ATH_CHECK( handle.isValid() );
615 trackCols.push_back(handle.cptr());
616 }
617
618 // output V0 vertices
619 SG::WriteHandle<xAOD::VertexContainer> V0OutputContainer;
620 if(m_vertexV0ContainerKey.key()=="" && m_v0VtxOutputKey.key()!="") {
621 V0OutputContainer = SG::WriteHandle<xAOD::VertexContainer>(m_v0VtxOutputKey, ctx);
622 ATH_CHECK( V0OutputContainer.record(std::make_unique<xAOD::VertexContainer>(), std::make_unique<xAOD::VertexAuxContainer>()) );
623 }
624
625 // Get Jpsi+X container
626 // Note: If the event does not contain a JX candidate, it is skipped and the V0 container will not be constructed if not previously available in StoreGate
627 SG::ReadHandle<xAOD::VertexContainer> jxContainer(m_vertexJXContainerKey, ctx);
628 ATH_CHECK( jxContainer.isValid() );
629 if(jxContainer->size()==0) return StatusCode::SUCCESS;
630
631 // Select the displaced tracks
632 std::vector<const xAOD::TrackParticle*> tracksDisplaced;
633 if(m_v0VtxOutputKey.key()!="" || m_disVDaug_num==3) {
634 for(const xAOD::TrackParticle* TP : *trackContainer.cptr()) {
635 // V0 track selection (https://gitlab.cern.ch/atlas/athena/-/blob/main/InnerDetector/InDetRecTools/InDetTrackSelectorTool/src/InDetConversionTrackSelectorTool.cxx)
636 if(m_v0TrkSelector->decision(*TP, primaryVertex)) {
637 uint8_t temp(0);
638 uint8_t nclus(0);
639 if(TP->summaryValue(temp, xAOD::numberOfPixelHits)) nclus += temp;
640 if(TP->summaryValue(temp, xAOD::numberOfSCTHits) ) nclus += temp;
641 if(!m_useTRT && nclus == 0) continue;
642
643 bool trk_cut = false;
644 if(nclus != 0) trk_cut = true;
645 if(nclus == 0 && TP->pt()>=m_ptTRT) trk_cut = true;
646 if(!trk_cut) continue;
647
648 // track is used if std::abs(d0/sig_d0) > d0_cut for PV
649 if(!d0Pass(ctx,TP,primaryVertex)) continue;
650
651 tracksDisplaced.push_back(TP);
652 }
653 }
654 }
655
656 SG::AuxElement::Accessor<std::string> mAcc_type("Type_V0Vtx");
657 SG::AuxElement::Accessor<int> mAcc_gfit("gamma_fit");
658 SG::AuxElement::Accessor<float> mAcc_gmass("gamma_mass");
659 SG::AuxElement::Accessor<float> mAcc_gchisq("gamma_chisq");
660 SG::AuxElement::Accessor<int> mAcc_gndof("gamma_ndof");
661
662 std::vector<std::pair<const xAOD::Vertex*,V0Enum> > selectedV0Candidates;
663
664 SG::ReadHandle<xAOD::VertexContainer> V0Container;
665 if(m_vertexV0ContainerKey.key() != "") {
666 V0Container = SG::ReadHandle<xAOD::VertexContainer>(m_vertexV0ContainerKey, ctx);
667 ATH_CHECK( V0Container.isValid() );
668
669 for(const xAOD::Vertex* vtx : *V0Container.cptr()) {
670 std::string type_V0Vtx;
671 if(mAcc_type.isAvailable(*vtx)) type_V0Vtx = mAcc_type(*vtx);
672
673 V0Enum opt(UNKNOWN); double massV0(0);
674 if(type_V0Vtx == "Lambda") {
675 opt = LAMBDA;
676 massV0 = m_V0Tools->invariantMass(vtx, m_massesV0_ppi);
677 if(massV0<m_LambdaMassLower || massV0>m_LambdaMassUpper) continue;
678 }
679 else if(type_V0Vtx == "Lambdabar") {
680 opt = LAMBDABAR;
681 massV0 = m_V0Tools->invariantMass(vtx, m_massesV0_pip);
682 if(massV0<m_LambdaMassLower || massV0>m_LambdaMassUpper) continue;
683 }
684 else if(type_V0Vtx == "Ks") {
685 opt = KS;
686 massV0 = m_V0Tools->invariantMass(vtx, m_massesV0_pipi);
687 if(massV0<m_KsMassLower || massV0>m_KsMassUpper) continue;
688 }
689
690 if(opt==UNKNOWN) continue;
691 if((opt==LAMBDA || opt==LAMBDABAR) && m_V0Hypothesis == "Ks") continue;
692 if(opt==KS && m_V0Hypothesis == "Lambda") continue;
693
694 int gamma_fit = mAcc_gfit.isAvailable(*vtx) ? mAcc_gfit(*vtx) : 0;
695 double gamma_mass = mAcc_gmass.isAvailable(*vtx) ? mAcc_gmass(*vtx) : -1;
696 double gamma_chisq = mAcc_gchisq.isAvailable(*vtx) ? mAcc_gchisq(*vtx) : 999999;
697 double gamma_ndof = mAcc_gndof.isAvailable(*vtx) ? mAcc_gndof(*vtx) : 0;
698 if(gamma_fit==1 && gamma_mass<m_minMass_gamma && gamma_chisq/gamma_ndof<m_chi2cut_gamma) continue;
699
700 selectedV0Candidates.push_back(std::pair<const xAOD::Vertex*,V0Enum>{vtx,opt});
701 }
702 }
703 else {
704 // fit V0 vertices
705 fitV0Container(ctx, V0OutputContainer.ptr(), tracksDisplaced, trackCols);
706
707 for(const xAOD::Vertex* vtx : *V0OutputContainer.cptr()) {
708 std::string type_V0Vtx;
709 if(mAcc_type.isAvailable(*vtx)) type_V0Vtx = mAcc_type(*vtx);
710
711 V0Enum opt(UNKNOWN); double massV0(0);
712 if(type_V0Vtx == "Lambda") {
713 opt = LAMBDA;
714 massV0 = m_V0Tools->invariantMass(vtx, m_massesV0_ppi);
715 if(massV0<m_LambdaMassLower || massV0>m_LambdaMassUpper) continue;
716 }
717 else if(type_V0Vtx == "Lambdabar") {
718 opt = LAMBDABAR;
719 massV0 = m_V0Tools->invariantMass(vtx, m_massesV0_pip);
720 if(massV0<m_LambdaMassLower || massV0>m_LambdaMassUpper) continue;
721 }
722 else if(type_V0Vtx == "Ks") {
723 opt = KS;
724 massV0 = m_V0Tools->invariantMass(vtx, m_massesV0_pipi);
725 if(massV0<m_KsMassLower || massV0>m_KsMassUpper) continue;
726 }
727
728 if(opt==UNKNOWN) continue;
729 if((opt==LAMBDA || opt==LAMBDABAR) && m_V0Hypothesis == "Ks") continue;
730 if(opt==KS && m_V0Hypothesis == "Lambda") continue;
731
732 int gamma_fit = mAcc_gfit.isAvailable(*vtx) ? mAcc_gfit(*vtx) : 0;
733 double gamma_mass = mAcc_gmass.isAvailable(*vtx) ? mAcc_gmass(*vtx) : -1;
734 double gamma_chisq = mAcc_gchisq.isAvailable(*vtx) ? mAcc_gchisq(*vtx) : 999999;
735 double gamma_ndof = mAcc_gndof.isAvailable(*vtx) ? mAcc_gndof(*vtx) : 0;
736 if(gamma_fit==1 && gamma_mass<m_minMass_gamma && gamma_chisq/gamma_ndof<m_chi2cut_gamma) continue;
737
738 selectedV0Candidates.push_back(std::pair<const xAOD::Vertex*,V0Enum>{vtx,opt});
739 }
740 }
741
742 // sort and chop the V0 candidates
743 std::sort( selectedV0Candidates.begin(), selectedV0Candidates.end(), [](std::pair<const xAOD::Vertex*,V0Enum>& a, std::pair<const xAOD::Vertex*,V0Enum>& b) { return a.first->chiSquared()/a.first->numberDoF() < b.first->chiSquared()/b.first->numberDoF(); } );
744 if(m_maxV0Candidates>0 && selectedV0Candidates.size()>m_maxV0Candidates) {
745 selectedV0Candidates.erase(selectedV0Candidates.begin()+m_maxV0Candidates, selectedV0Candidates.end());
746 }
747 if(selectedV0Candidates.size()==0) return StatusCode::SUCCESS;
748
749 std::vector<std::pair<Trk::VxCascadeInfo*, Trk::VxCascadeInfo*> > cascadeinfoContainer;
750 ATH_CHECK( performSearch(cascadeinfoContainer, selectedV0Candidates, tracksDisplaced, ctx) );
751
752 // sort and chop the main candidates
753 std::sort( cascadeinfoContainer.begin(), cascadeinfoContainer.end(), [](std::pair<Trk::VxCascadeInfo*, Trk::VxCascadeInfo*> a, std::pair<Trk::VxCascadeInfo*, Trk::VxCascadeInfo*> b) { return a.first->fitChi2()/a.first->nDoF() < b.first->fitChi2()/b.first->nDoF(); } );
754 if(m_maxMainVCandidates>0 && cascadeinfoContainer.size()>m_maxMainVCandidates) {
755 for(auto it=cascadeinfoContainer.begin()+m_maxMainVCandidates; it!=cascadeinfoContainer.end(); it++) {
756 if(it->first) delete it->first;
757 if(it->second) delete it->second;
758 }
759 cascadeinfoContainer.erase(cascadeinfoContainer.begin()+m_maxMainVCandidates, cascadeinfoContainer.end());
760 }
761 if(cascadeinfoContainer.size()==0) return StatusCode::SUCCESS;
762
763 SG::ReadHandle<xAOD::EventInfo> evt(m_eventInfo_key, ctx);
764 ATH_CHECK( evt.isValid() );
765 BPhysPVCascadeTools helper(&(*m_CascadeTools), evt.cptr());
766 helper.SetMinNTracksInPV(m_PV_minNTracks);
767
768 // Decorators for the cascade vertices
769 SG::AuxElement::Decorator<VertexLinkVector> CascadeLinksDecor("CascadeVertexLinks");
770 SG::AuxElement::Decorator<VertexLinkVector> PrecedingLinksDecor("PrecedingVertexLinks");
771 SG::AuxElement::Decorator<float> chi2_decor("ChiSquared");
772 SG::AuxElement::Decorator<int> ndof_decor("nDoF");
773 SG::AuxElement::Decorator<float> Pt_decor("Pt");
774 SG::AuxElement::Decorator<float> PtErr_decor("PtErr");
775
776 SG::AuxElement::Decorator<float> lxy_SV0_decor("lxy_SV0");
777 SG::AuxElement::Decorator<float> lxyErr_SV0_decor("lxyErr_SV0");
778 SG::AuxElement::Decorator<float> a0xy_SV0_decor("a0xy_SV0");
779 SG::AuxElement::Decorator<float> a0xyErr_SV0_decor("a0xyErr_SV0");
780 SG::AuxElement::Decorator<float> a0z_SV0_decor("a0z_SV0");
781 SG::AuxElement::Decorator<float> a0zErr_SV0_decor("a0zErr_SV0");
782
783 SG::AuxElement::Decorator<float> lxy_SV1_decor("lxy_SV1");
784 SG::AuxElement::Decorator<float> lxyErr_SV1_decor("lxyErr_SV1");
785 SG::AuxElement::Decorator<float> a0xy_SV1_decor("a0xy_SV1");
786 SG::AuxElement::Decorator<float> a0xyErr_SV1_decor("a0xyErr_SV1");
787 SG::AuxElement::Decorator<float> a0z_SV1_decor("a0z_SV1");
788 SG::AuxElement::Decorator<float> a0zErr_SV1_decor("a0zErr_SV1");
789
790 SG::AuxElement::Decorator<float> lxy_SV2_decor("lxy_SV2");
791 SG::AuxElement::Decorator<float> lxyErr_SV2_decor("lxyErr_SV2");
792 SG::AuxElement::Decorator<float> a0xy_SV2_decor("a0xy_SV2");
793 SG::AuxElement::Decorator<float> a0xyErr_SV2_decor("a0xyErr_SV2");
794 SG::AuxElement::Decorator<float> a0z_SV2_decor("a0z_SV2");
795 SG::AuxElement::Decorator<float> a0zErr_SV2_decor("a0zErr_SV2");
796
797 SG::AuxElement::Decorator<float> chi2_V2_decor("ChiSquared_V2");
798 SG::AuxElement::Decorator<int> ndof_V2_decor("nDoF_V2");
799
800 for(auto cascade_info_pair : cascadeinfoContainer) {
801 if(cascade_info_pair.first==nullptr) {
802 ATH_MSG_ERROR("CascadeInfo is null");
803 continue;
804 }
805
806 const std::vector<xAOD::Vertex*> &cascadeVertices = cascade_info_pair.first->vertices();
807 if(cascadeVertices.size() != topoN) ATH_MSG_ERROR("Incorrect number of vertices");
808 for(size_t i=0; i<topoN; i++) {
809 if(cascadeVertices[i]==nullptr) ATH_MSG_ERROR("Error null vertex");
810 }
811
812 cascade_info_pair.first->setSVOwnership(false); // Prevent Container from deleting vertices
813 const auto mainVertex = cascadeVertices[topoN-1]; // this is the mother vertex
814 const std::vector< std::vector<TLorentzVector> > &moms = cascade_info_pair.first->getParticleMoms();
815
816 // Identify the input JX
817 int ijx = m_JXSubVtx ? topoN-2 : topoN-1;
818 if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0) {
819 if(m_JXV0SubVtx) ijx = 1;
820 else ijx = topoN-1;
821 }
822 else if(m_extraTrk1MassHypo>0 && m_disVDaug_num==2) {
823 if(m_JXV0SubVtx || m_JXSubVtx) ijx = 1;
824 else ijx = topoN-1;
825 }
826 const xAOD::Vertex* jxVtx(nullptr);
827 if(m_jxDaug_num==4) jxVtx = FindVertex<4>(jxContainer.ptr(), cascadeVertices[ijx]);
828 else if(m_jxDaug_num==3) jxVtx = FindVertex<3>(jxContainer.ptr(), cascadeVertices[ijx]);
829 else jxVtx = FindVertex<2>(jxContainer.ptr(), cascadeVertices[ijx]);
830
831 xAOD::BPhysHypoHelper vtx(m_hypoName, mainVertex);
832
833 // Get refitted track momenta from all vertices, charged tracks only
834 BPhysPVCascadeTools::SetVectorInfo(vtx, cascade_info_pair.first);
835 vtx.setPass(true);
836
837 //
838 // Decorate main vertex
839 //
840 // mass, mass error
841 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/VxCascadeInfo.h
842 BPHYS_CHECK( vtx.setMass(m_CascadeTools->invariantMass(moms[topoN-1])) );
843 BPHYS_CHECK( vtx.setMassErr(m_CascadeTools->invariantMassError(moms[topoN-1],cascade_info_pair.first->getCovariance()[topoN-1])) );
844 // pt and pT error (the default pt of mainVertex is != the pt of the full cascade fit!)
845 Pt_decor(*mainVertex) = m_CascadeTools->pT(moms[topoN-1]);
846 PtErr_decor(*mainVertex) = m_CascadeTools->pTError(moms[topoN-1],cascade_info_pair.first->getCovariance()[topoN-1]);
847 // chi2 and ndof (the default chi2 of mainVertex is != the chi2 of the full cascade fit!)
848 chi2_decor(*mainVertex) = cascade_info_pair.first->fitChi2();
849 ndof_decor(*mainVertex) = cascade_info_pair.first->nDoF();
850
851 if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0) { // special cases
852 if(m_JXV0SubVtx) {
853 lxy_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
854 lxyErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->lxyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
855 a0z_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0z(moms[0],cascadeVertices[0],cascadeVertices[1]);
856 a0zErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0zError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
857 a0xy_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0xy(moms[0],cascadeVertices[0],cascadeVertices[1]);
858 a0xyErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0xyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
859 lxy_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->lxy(moms[1],cascadeVertices[1],mainVertex);
860 lxyErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->lxyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
861 a0z_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0z(moms[1],cascadeVertices[1],mainVertex);
862 a0zErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0zError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
863 a0xy_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0xy(moms[1],cascadeVertices[1],mainVertex);
864 a0xyErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0xyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
865 lxy_SV2_decor(*cascadeVertices[2]) = m_CascadeTools->lxy(moms[2],cascadeVertices[2],mainVertex);
866 lxyErr_SV2_decor(*cascadeVertices[2]) = m_CascadeTools->lxyError(moms[2],cascade_info_pair.first->getCovariance()[2],cascadeVertices[2],mainVertex);
867 a0z_SV2_decor(*cascadeVertices[2]) = m_CascadeTools->a0z(moms[2],cascadeVertices[2],mainVertex);
868 a0zErr_SV2_decor(*cascadeVertices[2]) = m_CascadeTools->a0zError(moms[2],cascade_info_pair.first->getCovariance()[2],cascadeVertices[2],mainVertex);
869 a0xy_SV2_decor(*cascadeVertices[2]) = m_CascadeTools->a0xy(moms[2],cascadeVertices[2],mainVertex);
870 a0xyErr_SV2_decor(*cascadeVertices[2]) = m_CascadeTools->a0xyError(moms[2],cascade_info_pair.first->getCovariance()[2],cascadeVertices[2],mainVertex);
871 }
872 else {
873 lxy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],mainVertex);
874 lxyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],mainVertex);
875 a0z_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0z(moms[0],cascadeVertices[0],mainVertex);
876 a0zErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0zError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],mainVertex);
877 a0xy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xy(moms[0],cascadeVertices[0],mainVertex);
878 a0xyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],mainVertex);
879 lxy_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->lxy(moms[1],cascadeVertices[1],mainVertex);
880 lxyErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->lxyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
881 a0z_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0z(moms[1],cascadeVertices[1],mainVertex);
882 a0zErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0zError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
883 a0xy_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0xy(moms[1],cascadeVertices[1],mainVertex);
884 a0xyErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0xyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
885 }
886 }
887 else if(m_extraTrk1MassHypo>0 && m_disVDaug_num==2) { // special cases
888 if(m_JXV0SubVtx) {
889 lxy_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
890 lxyErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->lxyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
891 a0z_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0z(moms[0],cascadeVertices[0],cascadeVertices[1]);
892 a0zErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0zError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
893 a0xy_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0xy(moms[0],cascadeVertices[0],cascadeVertices[1]);
894 a0xyErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0xyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
895 lxy_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->lxy(moms[1],cascadeVertices[1],mainVertex);
896 lxyErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->lxyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
897 a0z_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0z(moms[1],cascadeVertices[1],mainVertex);
898 a0zErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0zError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
899 a0xy_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0xy(moms[1],cascadeVertices[1],mainVertex);
900 a0xyErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0xyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
901 }
902 else {
903 lxy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],mainVertex);
904 lxyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],mainVertex);
905 a0z_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0z(moms[0],cascadeVertices[0],mainVertex);
906 a0zErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0zError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],mainVertex);
907 a0xy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xy(moms[0],cascadeVertices[0],mainVertex);
908 a0xyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],mainVertex);
909 if(m_JXSubVtx) {
910 lxy_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->lxy(moms[1],cascadeVertices[1],mainVertex);
911 lxyErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->lxyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
912 a0z_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0z(moms[1],cascadeVertices[1],mainVertex);
913 a0zErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0zError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
914 a0xy_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0xy(moms[1],cascadeVertices[1],mainVertex);
915 a0xyErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0xyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
916 }
917 }
918 }
919 else { // other normal cases
920 if(m_disVDaug_num==2) {
921 lxy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],mainVertex);
922 lxyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],mainVertex);
923 a0z_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0z(moms[0],cascadeVertices[0],mainVertex);
924 a0zErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0zError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],mainVertex);
925 a0xy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xy(moms[0],cascadeVertices[0],mainVertex);
926 a0xyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],mainVertex);
927 }
928 else {
929 lxy_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
930 lxyErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->lxyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
931 a0z_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0z(moms[0],cascadeVertices[0],cascadeVertices[1]);
932 a0zErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0zError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
933 a0xy_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0xy(moms[0],cascadeVertices[0],cascadeVertices[1]);
934 a0xyErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0xyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
935 lxy_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->lxy(moms[1],cascadeVertices[1],mainVertex);
936 lxyErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->lxyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
937 a0xy_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0z(moms[1],cascadeVertices[1],mainVertex);
938 a0xyErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0zError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
939 a0z_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0xy(moms[1],cascadeVertices[1],mainVertex);
940 a0zErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0xyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
941 }
942
943 if(m_JXSubVtx) {
944 lxy_SV2_decor(*cascadeVertices[ijx]) = m_CascadeTools->lxy(moms[ijx],cascadeVertices[ijx],mainVertex);
945 lxyErr_SV2_decor(*cascadeVertices[ijx]) = m_CascadeTools->lxyError(moms[ijx],cascade_info_pair.first->getCovariance()[ijx],cascadeVertices[ijx],mainVertex);
946 a0z_SV2_decor(*cascadeVertices[ijx]) = m_CascadeTools->a0z(moms[ijx],cascadeVertices[ijx],mainVertex);
947 a0zErr_SV2_decor(*cascadeVertices[ijx]) = m_CascadeTools->a0zError(moms[ijx],cascade_info_pair.first->getCovariance()[ijx],cascadeVertices[ijx],mainVertex);
948 a0xy_SV2_decor(*cascadeVertices[ijx]) = m_CascadeTools->a0xy(moms[ijx],cascadeVertices[ijx],mainVertex);
949 a0xyErr_SV2_decor(*cascadeVertices[ijx]) = m_CascadeTools->a0xyError(moms[ijx],cascade_info_pair.first->getCovariance()[ijx],cascadeVertices[ijx],mainVertex);
950 }
951 }
952
953 chi2_V2_decor(*cascadeVertices[ijx]) = m_V0Tools->chisq(jxVtx);
954 ndof_V2_decor(*cascadeVertices[ijx]) = m_V0Tools->ndof(jxVtx);
955
956 if(m_cascadeFitWithPV==0) {
957 ATH_CHECK(helper.FillCandwithRefittedVertices(m_refitPV, defaultPVContainer.cptr(), m_refitPV ? refPvContainer.ptr() : 0, &(*m_pvRefitter), m_PV_max, m_DoVertexType, cascade_info_pair.first, topoN-1, m_massMainV, vtx));
958 }
959
960 for(size_t i=0; i<topoN; i++) {
961 VtxWriteHandles[i].ptr()->push_back(cascadeVertices[i]);
962 }
963
964 // Set links to cascade vertices
965 VertexLinkVector cascadeVertexLinks;
966 VertexLink vertexLink1;
967 vertexLink1.setElement(cascadeVertices[0]);
968 vertexLink1.setStorableObject(*VtxWriteHandles[0].ptr());
969 if( vertexLink1.isValid() ) cascadeVertexLinks.push_back( vertexLink1 );
970 if(topoN>=3) {
971 VertexLink vertexLink2;
972 vertexLink2.setElement(cascadeVertices[1]);
973 vertexLink2.setStorableObject(*VtxWriteHandles[1].ptr());
974 if( vertexLink2.isValid() ) cascadeVertexLinks.push_back( vertexLink2 );
975 }
976 if(topoN==4) {
977 VertexLink vertexLink3;
978 vertexLink3.setElement(cascadeVertices[2]);
979 vertexLink3.setStorableObject(*VtxWriteHandles[2].ptr());
980 if( vertexLink3.isValid() ) cascadeVertexLinks.push_back( vertexLink3 );
981 }
982 CascadeLinksDecor(*mainVertex) = cascadeVertexLinks;
983
984 // process the posterior cascade with mass constraint for the main vertex
985 if(cascade_info_pair.second) {
986 const std::vector<xAOD::Vertex*> &cascadeVertices_mvc = cascade_info_pair.second->vertices();
987 if(cascadeVertices_mvc.size() != topoN) ATH_MSG_ERROR("Incorrect number of vertices (mvc)");
988 for(size_t i=0; i<topoN; i++) {
989 if(cascadeVertices_mvc[i]==nullptr) ATH_MSG_ERROR("Error null vertex (mvc)");
990 }
991 cascade_info_pair.second->setSVOwnership(false);
992 const auto mainVertex_mvc = cascadeVertices_mvc[topoN-1];
993 const std::vector< std::vector<TLorentzVector> > &moms_mvc = cascade_info_pair.second->getParticleMoms();
994 // Identify the input JX
995 int ijx_mvc = m_JXSubVtx ? topoN-2 : topoN-1;
996 if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0) {
997 if(m_JXV0SubVtx) ijx_mvc = 1;
998 else ijx_mvc = topoN-1;
999 }
1000 else if(m_extraTrk1MassHypo>0 && m_disVDaug_num==2) {
1001 if(m_JXV0SubVtx || m_JXSubVtx) ijx_mvc = 1;
1002 else ijx_mvc = topoN-1;
1003 }
1004 const xAOD::Vertex* jxVtx_mvc(nullptr);
1005 if(m_jxDaug_num==4) jxVtx_mvc = FindVertex<4>(jxContainer.ptr(), cascadeVertices_mvc[ijx_mvc]);
1006 else if(m_jxDaug_num==3) jxVtx_mvc = FindVertex<3>(jxContainer.ptr(), cascadeVertices_mvc[ijx_mvc]);
1007 else jxVtx_mvc = FindVertex<2>(jxContainer.ptr(), cascadeVertices_mvc[ijx_mvc]);
1008
1009 xAOD::BPhysHypoHelper vtx_mvc(m_hypoName, mainVertex_mvc);
1010 BPhysPVCascadeTools::SetVectorInfo(vtx_mvc, cascade_info_pair.second);
1011 vtx_mvc.setPass(true);
1012
1013 BPHYS_CHECK( vtx_mvc.setMass(m_CascadeTools->invariantMass(moms_mvc[topoN-1])) );
1014 BPHYS_CHECK( vtx_mvc.setMassErr(m_CascadeTools->invariantMassError(moms_mvc[topoN-1],cascade_info_pair.second->getCovariance()[topoN-1])) );
1015 Pt_decor(*mainVertex_mvc) = m_CascadeTools->pT(moms_mvc[topoN-1]);
1016 PtErr_decor(*mainVertex_mvc) = m_CascadeTools->pTError(moms_mvc[topoN-1],cascade_info_pair.second->getCovariance()[topoN-1]);
1017 chi2_decor(*mainVertex_mvc) = cascade_info_pair.second->fitChi2();
1018 ndof_decor(*mainVertex_mvc) = cascade_info_pair.second->nDoF();
1019
1020 if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0) { // special cases
1021 if(m_JXV0SubVtx) {
1022 lxy_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxy(moms_mvc[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1023 lxyErr_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1024 a0z_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0z(moms_mvc[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1025 a0zErr_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0zError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1026 a0xy_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xy(moms_mvc[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1027 a0xyErr_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1028 lxy_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1029 lxyErr_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1030 a0z_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0z(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1031 a0zErr_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0zError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1032 a0xy_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1033 a0xyErr_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1034 lxy_SV2_decor(*cascadeVertices_mvc[2]) = m_CascadeTools->lxy(moms_mvc[2],cascadeVertices_mvc[2],mainVertex_mvc);
1035 lxyErr_SV2_decor(*cascadeVertices_mvc[2]) = m_CascadeTools->lxyError(moms_mvc[2],cascade_info_pair.second->getCovariance()[2],cascadeVertices_mvc[2],mainVertex_mvc);
1036 a0z_SV2_decor(*cascadeVertices_mvc[2]) = m_CascadeTools->a0z(moms_mvc[2],cascadeVertices_mvc[2],mainVertex_mvc);
1037 a0zErr_SV2_decor(*cascadeVertices_mvc[2]) = m_CascadeTools->a0zError(moms_mvc[2],cascade_info_pair.second->getCovariance()[2],cascadeVertices_mvc[2],mainVertex_mvc);
1038 a0xy_SV2_decor(*cascadeVertices_mvc[2]) = m_CascadeTools->a0xy(moms_mvc[2],cascadeVertices_mvc[2],mainVertex_mvc);
1039 a0xyErr_SV2_decor(*cascadeVertices_mvc[2]) = m_CascadeTools->a0xyError(moms_mvc[2],cascade_info_pair.second->getCovariance()[2],cascadeVertices_mvc[2],mainVertex_mvc);
1040 }
1041 else {
1042 lxy_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxy(moms_mvc[0],cascadeVertices_mvc[0],mainVertex_mvc);
1043 lxyErr_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],mainVertex_mvc);
1044 a0z_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0z(moms_mvc[0],cascadeVertices_mvc[0],mainVertex_mvc);
1045 a0zErr_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0zError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],mainVertex_mvc);
1046 a0xy_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xy(moms_mvc[0],cascadeVertices_mvc[0],mainVertex_mvc);
1047 a0xyErr_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],mainVertex_mvc);
1048 lxy_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1049 lxyErr_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1050 a0z_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0z(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1051 a0zErr_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0zError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1052 a0xy_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1053 a0xyErr_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1054 }
1055 }
1056 else if(m_extraTrk1MassHypo>0 && m_disVDaug_num==2) { // special cases
1057 if(m_JXV0SubVtx) {
1058 lxy_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxy(moms_mvc[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1059 lxyErr_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1060 a0z_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0z(moms_mvc[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1061 a0zErr_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0zError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1062 a0xy_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xy(moms_mvc[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1063 a0xyErr_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1064 lxy_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1065 lxyErr_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1066 a0z_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0z(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1067 a0zErr_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0zError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1068 a0xy_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1069 a0xyErr_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1070 }
1071 else {
1072 lxy_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxy(moms_mvc[0],cascadeVertices_mvc[0],mainVertex_mvc);
1073 lxyErr_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],mainVertex_mvc);
1074 a0z_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0z(moms_mvc[0],cascadeVertices_mvc[0],mainVertex_mvc);
1075 a0zErr_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0zError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],mainVertex_mvc);
1076 a0xy_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xy(moms_mvc[0],cascadeVertices_mvc[0],mainVertex_mvc);
1077 a0xyErr_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],mainVertex_mvc);
1078 if(m_JXSubVtx) {
1079 lxy_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1080 lxyErr_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1081 a0z_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0z(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1082 a0zErr_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0zError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1083 a0xy_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1084 a0xyErr_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1085 }
1086 }
1087 }
1088 else { // other normal cases
1089 if(m_disVDaug_num==2) {
1090 lxy_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxy(moms_mvc[0],cascadeVertices_mvc[0],mainVertex_mvc);
1091 lxyErr_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],mainVertex_mvc);
1092 a0z_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0z(moms_mvc[0],cascadeVertices_mvc[0],mainVertex_mvc);
1093 a0zErr_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0zError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],mainVertex_mvc);
1094 a0xy_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xy(moms_mvc[0],cascadeVertices_mvc[0],mainVertex_mvc);
1095 a0xyErr_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],mainVertex_mvc);
1096 }
1097 else {
1098 lxy_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxy(moms_mvc[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1099 lxyErr_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1100 a0z_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0z(moms_mvc[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1101 a0zErr_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0zError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1102 a0xy_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xy(moms_mvc[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1103 a0xyErr_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1104 lxy_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1105 lxyErr_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1106 a0xy_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0z(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1107 a0xyErr_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0zError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1108 a0z_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1109 a0zErr_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1110 }
1111 if(m_JXSubVtx) {
1112 lxy_SV2_decor(*cascadeVertices_mvc[ijx_mvc]) = m_CascadeTools->lxy(moms_mvc[ijx_mvc],cascadeVertices_mvc[ijx_mvc],mainVertex_mvc);
1113 lxyErr_SV2_decor(*cascadeVertices_mvc[ijx_mvc]) = m_CascadeTools->lxyError(moms_mvc[ijx_mvc],cascade_info_pair.second->getCovariance()[ijx_mvc],cascadeVertices_mvc[ijx_mvc],mainVertex_mvc);
1114 a0z_SV2_decor(*cascadeVertices_mvc[ijx_mvc]) = m_CascadeTools->a0z(moms_mvc[ijx_mvc],cascadeVertices_mvc[ijx_mvc],mainVertex_mvc);
1115 a0zErr_SV2_decor(*cascadeVertices_mvc[ijx_mvc]) = m_CascadeTools->a0zError(moms_mvc[ijx_mvc],cascade_info_pair.second->getCovariance()[ijx_mvc],cascadeVertices_mvc[ijx_mvc],mainVertex_mvc);
1116 a0xy_SV2_decor(*cascadeVertices_mvc[ijx_mvc]) = m_CascadeTools->a0xy(moms_mvc[ijx_mvc],cascadeVertices_mvc[ijx_mvc],mainVertex_mvc);
1117 a0xyErr_SV2_decor(*cascadeVertices_mvc[ijx_mvc]) = m_CascadeTools->a0xyError(moms_mvc[ijx_mvc],cascade_info_pair.second->getCovariance()[ijx_mvc],cascadeVertices_mvc[ijx_mvc],mainVertex_mvc);
1118 }
1119 }
1120 chi2_V2_decor(*cascadeVertices_mvc[ijx_mvc]) = m_V0Tools->chisq(jxVtx_mvc);
1121 ndof_V2_decor(*cascadeVertices_mvc[ijx_mvc]) = m_V0Tools->ndof(jxVtx_mvc);
1122
1123 if(m_cascadeFitWithPV==0) {
1124 ATH_CHECK(helper.FillCandwithRefittedVertices(m_refitPV, defaultPVContainer.cptr(), m_refitPV ? refPvContainer.ptr() : 0, &(*m_pvRefitter), m_PV_max, m_DoVertexType, cascade_info_pair.second, topoN-1, m_massMainV, vtx_mvc));
1125 }
1126
1127 for(size_t i=0; i<topoN; i++) {
1128 VtxWriteHandles_mvc[i].ptr()->push_back(cascadeVertices_mvc[i]);
1129 }
1130
1131 VertexLinkVector cascadeVertexLinks_mvc;
1132 VertexLink vertexLink1_mvc;
1133 vertexLink1_mvc.setElement(cascadeVertices_mvc[0]);
1134 vertexLink1_mvc.setStorableObject(*VtxWriteHandles_mvc[0].ptr());
1135 if( vertexLink1_mvc.isValid() ) cascadeVertexLinks_mvc.push_back( vertexLink1_mvc );
1136 if(topoN>=3) {
1137 VertexLink vertexLink2_mvc;
1138 vertexLink2_mvc.setElement(cascadeVertices_mvc[1]);
1139 vertexLink2_mvc.setStorableObject(*VtxWriteHandles_mvc[1].ptr());
1140 if( vertexLink2_mvc.isValid() ) cascadeVertexLinks_mvc.push_back( vertexLink2_mvc );
1141 }
1142 if(topoN==4) {
1143 VertexLink vertexLink3_mvc;
1144 vertexLink3_mvc.setElement(cascadeVertices_mvc[2]);
1145 vertexLink3_mvc.setStorableObject(*VtxWriteHandles_mvc[2].ptr());
1146 if( vertexLink3_mvc.isValid() ) cascadeVertexLinks_mvc.push_back( vertexLink3_mvc );
1147 }
1148 CascadeLinksDecor(*mainVertex_mvc) = cascadeVertexLinks_mvc;
1149
1150 VertexLinkVector precedingVertexLinks_mvc;
1151 VertexLink vertexLink_mvc;
1152 vertexLink_mvc.setElement(cascadeVertices[topoN-1]);
1153 vertexLink_mvc.setStorableObject(*VtxWriteHandles[topoN-1].ptr());
1154 if( vertexLink_mvc.isValid() ) precedingVertexLinks_mvc.push_back( vertexLink_mvc );
1155 PrecedingLinksDecor(*mainVertex_mvc) = precedingVertexLinks_mvc;
1156 } // cascade_info_pair.second!=0
1157 } // loop over cascadeinfoContainer
1158
1159 // Deleting cascadeinfo since this won't be stored.
1160 for (auto cascade_info_pair : cascadeinfoContainer) {
1161 if(cascade_info_pair.first) delete cascade_info_pair.first;
1162 if(cascade_info_pair.second) delete cascade_info_pair.second;
1163 }
1164
1165 return StatusCode::SUCCESS;
1166 }
1167
1168 bool JpsiXPlusDisplaced::d0Pass(const EventContext& ctx, const xAOD::TrackParticle* track, const xAOD::Vertex* PV) const {
1169 bool pass = false;
1170 std::unique_ptr<Trk::Perigee> per = m_trackToVertexTool->perigeeAtVertex(ctx, *track, PV->position());
1171 if(!per) return pass;
1172 double d0 = per->parameters()[Trk::d0];
1173 double sig_d0 = sqrt((*per->covariance())(0,0));
1174 if(std::abs(d0/sig_d0) > m_d0_cut) pass = true;
1175 return pass;
1176 }
1177
1178 JpsiXPlusDisplaced::XiCandidate JpsiXPlusDisplaced::getXiCandidate(const EventContext& ctx, const xAOD::Vertex* V0vtx, const V0Enum V0, const xAOD::TrackParticle* track3) const {
1179 XiCandidate disVtx;
1180
1181 std::vector<const xAOD::TrackParticle*> tracksV0;
1182 tracksV0.reserve(V0vtx->nTrackParticles());
1183 for(size_t i=0; i<V0vtx->nTrackParticles(); i++) tracksV0.push_back(V0vtx->trackParticle(i));
1184 std::vector<double> massesV0;
1185 if(V0==LAMBDA) massesV0 = m_massesV0_ppi;
1186 else if(V0==LAMBDABAR) massesV0 = m_massesV0_pip;
1187 else if(V0==KS) massesV0 = m_massesV0_pipi;
1188
1189 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
1190 int robustness = 0;
1191 m_iVertexFitter->setRobustness(robustness, *state);
1192 std::vector<Trk::VertexID> vrtList;
1193 // V0 vertex
1194 Trk::VertexID vID = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1195 vrtList.push_back(vID);
1196 // Mother vertex
1197 std::vector<const xAOD::TrackParticle*> tracksDis3{track3};
1198 std::vector<double> massesDis3{m_disVDaug3MassHypo};
1199 m_iVertexFitter->nextVertex(tracksDis3,massesDis3,vrtList,*state);
1200 // Do the work
1201 std::unique_ptr<Trk::VxCascadeInfo> cascade_info = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state) );
1202 if(cascade_info) {
1203 cascade_info->setSVOwnership(true);
1204 double chi2NDF = cascade_info->fitChi2()/cascade_info->nDoF();
1205 if(m_chi2cut_DisV<=0 || chi2NDF < m_chi2cut_DisV) {
1206 const std::vector<std::vector<TLorentzVector> > &moms = cascade_info->getParticleMoms();
1207 const std::vector<xAOD::Vertex*> &cascadeVertices = cascade_info->vertices();
1208 double lxy_SV1_sub = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
1209 if(lxy_SV1_sub > 0.2) {
1210 TLorentzVector totalMom;
1211 for(size_t it=0; it<moms[1].size(); it++) totalMom += moms[1][it];
1212 double disV_mass = totalMom.M();
1213 if(m_useImprovedMass) {
1214 if((V0==LAMBDA || V0==LAMBDABAR) && m_massLd>0) disV_mass += - moms[1][1].M() + m_massLd;
1215 else if(V0==KS && m_massKs>0) disV_mass += - moms[1][1].M() + m_massKs;
1216 }
1217 if(disV_mass>m_DisplacedMassLower && disV_mass<m_DisplacedMassUpper) {
1218 disVtx.track = track3; disVtx.V0vtx = V0vtx; disVtx.V0type = V0;
1219 disVtx.chi2NDF = chi2NDF; disVtx.p4_V0track1 = moms[0][0];
1220 disVtx.p4_V0track2 = moms[0][1]; disVtx.p4_disVtrack = moms[1][0];
1221 }
1222 }
1223 }
1224 }
1225 return disVtx;
1226 }
1227
1228 std::unique_ptr<xAOD::Vertex> JpsiXPlusDisplaced::fitTracks(const EventContext& ctx, const xAOD::TrackParticle* track1, const xAOD::TrackParticle* track2, const xAOD::TrackParticle* track3) const {
1229 // Starting point
1230 const Trk::Perigee& aPerigee1 = track1->perigeeParameters();
1231 const Trk::Perigee& aPerigee2 = track2->perigeeParameters();
1232 int sflag(0), errorcode(0);
1233 Amg::Vector3D startingPoint = m_vertexEstimator->getCirclesIntersectionPoint(&aPerigee1,&aPerigee2,sflag,errorcode);
1234 if(errorcode) startingPoint(0) = startingPoint(1) = startingPoint(2) = 0.0;
1235 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
1236 // do the fit
1237 if(track3) {
1238 return m_iVertexFitter->fit(std::vector<const xAOD::TrackParticle*>{track1,track2,track3}, startingPoint, *state);
1239 }
1240 else {
1241 return m_iVertexFitter->fit(std::vector<const xAOD::TrackParticle*>{track1,track2}, startingPoint, *state);
1242 }
1243 }
1244
1245 JpsiXPlusDisplaced::MesonCandidate JpsiXPlusDisplaced::getDpmCandidate(const EventContext& ctx, const xAOD::Vertex* JXvtx, const xAOD::TrackParticle* extraTrk1, const xAOD::TrackParticle* extraTrk2, const xAOD::TrackParticle* extraTrk3) const {
1246 MesonCandidate Dpm;
1247 // Check overlap
1248 std::vector<const xAOD::TrackParticle*> tracksJX;
1249 tracksJX.reserve(JXvtx->nTrackParticles());
1250 for(size_t i=0; i<JXvtx->nTrackParticles(); i++) tracksJX.push_back(JXvtx->trackParticle(i));
1251
1252 TLorentzVector tmp1, tmp2, tmp3;
1253 tmp1.SetPtEtaPhiM(extraTrk1->pt(),extraTrk1->eta(),extraTrk1->phi(),m_extraTrk1MassHypo);
1254 tmp2.SetPtEtaPhiM(extraTrk2->pt(),extraTrk2->eta(),extraTrk2->phi(),m_extraTrk2MassHypo);
1255 tmp3.SetPtEtaPhiM(extraTrk3->pt(),extraTrk3->eta(),extraTrk3->phi(),m_extraTrk3MassHypo);
1256 if((tmp1+tmp2+tmp3).M() < m_DpmMassLower || (tmp1+tmp2+tmp3).M() > m_DpmMassUpper) return Dpm;
1257
1258 std::unique_ptr<xAOD::Vertex> vtx = fitTracks(ctx, extraTrk1, extraTrk2, extraTrk3);
1259 if(vtx) {
1260 double chi2NDF = vtx->chiSquared()/vtx->numberDoF();
1261 if(m_chi2cut_Dpm<=0.0 || chi2NDF < m_chi2cut_Dpm) {
1262 double lxyDpm = m_V0Tools->lxy(vtx.get(),JXvtx);
1263 if(lxyDpm>m_lxyDpm_cut) {
1264 Dpm.extraTrack1 = extraTrk1; Dpm.extraTrack2 = extraTrk2; Dpm.extraTrack3 = extraTrk3;
1265 Dpm.chi2NDF = chi2NDF;
1266 TVector3 tot_pt; TVector3 tmp;
1267 for(size_t i=0; i<vtx->vxTrackAtVertex().size(); ++i) {
1268 const Trk::TrackParameters* aPerigee = vtx->vxTrackAtVertex()[i].perigeeAtVertex();
1269 if(aPerigee) {
1270 tmp.SetXYZ(aPerigee->momentum()[Trk::px],aPerigee->momentum()[Trk::py],aPerigee->momentum()[Trk::pz]);
1271 tot_pt += tmp;
1272 }
1273 }
1274 Dpm.pt = tot_pt.Pt();
1275 }
1276 }
1277 }
1278 return Dpm;
1279 }
1280
1281 JpsiXPlusDisplaced::MesonCandidate JpsiXPlusDisplaced::getD0Candidate(const EventContext& ctx, const xAOD::Vertex* JXvtx, const xAOD::TrackParticle* extraTrk1, const xAOD::TrackParticle* extraTrk2) const {
1282 MesonCandidate D0;
1283
1284 TLorentzVector tmp1, tmp2;
1285 tmp1.SetPtEtaPhiM(extraTrk1->pt(),extraTrk1->eta(),extraTrk1->phi(),m_extraTrk1MassHypo);
1286 tmp2.SetPtEtaPhiM(extraTrk2->pt(),extraTrk2->eta(),extraTrk2->phi(),m_extraTrk2MassHypo);
1287 if((tmp1+tmp2).M() < m_D0MassLower || (tmp1+tmp2).M() > m_D0MassUpper) return D0;
1288
1289 std::unique_ptr<xAOD::Vertex> vtx = fitTracks(ctx, extraTrk1, extraTrk2);
1290 if(vtx) {
1291 double chi2NDF = vtx->chiSquared()/vtx->numberDoF();
1292 if(m_chi2cut_D0<=0.0 || chi2NDF < m_chi2cut_D0) {
1293 double lxyD0 = m_V0Tools->lxy(vtx.get(),JXvtx);
1294 if(lxyD0>m_lxyD0_cut) {
1295 D0.extraTrack1 = extraTrk1; D0.extraTrack2 = extraTrk2;
1296 D0.chi2NDF = chi2NDF;
1297 TVector3 tot_pt; TVector3 tmp;
1298 for(size_t i=0; i<vtx->vxTrackAtVertex().size(); ++i) {
1299 const Trk::TrackParameters* aPerigee = vtx->vxTrackAtVertex()[i].perigeeAtVertex();
1300 if(aPerigee) {
1301 tmp.SetXYZ(aPerigee->momentum()[Trk::px],aPerigee->momentum()[Trk::py],aPerigee->momentum()[Trk::pz]);
1302 tot_pt += tmp;
1303 }
1304 }
1305 D0.pt = tot_pt.Pt();
1306 }
1307 }
1308 }
1309 return D0;
1310 }
1311
1312 std::vector<std::pair<Trk::VxCascadeInfo*,Trk::VxCascadeInfo*> > JpsiXPlusDisplaced::fitMainVtx(const EventContext& ctx, const xAOD::Vertex* JXvtx, const std::vector<double>& massesJX, const xAOD::Vertex* V0vtx, const V0Enum V0, const xAOD::TrackParticleContainer* trackContainer, const std::vector<const xAOD::TrackParticleContainer*>& trackCols, const xAOD::VertexContainer* defaultPVContainer, const xAOD::VertexContainer* pvContainer) const {
1313 std::vector<std::pair<Trk::VxCascadeInfo*,Trk::VxCascadeInfo*> > result;
1314
1315 std::vector<const xAOD::TrackParticle*> tracksJX;
1316 tracksJX.reserve(JXvtx->nTrackParticles());
1317 for(size_t i=0; i<JXvtx->nTrackParticles(); i++) tracksJX.push_back(JXvtx->trackParticle(i));
1318 if (tracksJX.size() != massesJX.size()) {
1319 ATH_MSG_ERROR("Problems with JX input: number of tracks or track mass inputs is not correct!");
1320 return result;
1321 }
1322 // Check identical tracks in input
1323 if(std::find(tracksJX.cbegin(), tracksJX.cend(), V0vtx->trackParticle(0)) != tracksJX.cend()) return result;
1324 if(std::find(tracksJX.cbegin(), tracksJX.cend(), V0vtx->trackParticle(1)) != tracksJX.cend()) return result;
1325 std::vector<const xAOD::TrackParticle*> tracksV0;
1326 tracksV0.reserve(V0vtx->nTrackParticles());
1327 for(size_t j=0; j<V0vtx->nTrackParticles(); j++) tracksV0.push_back(V0vtx->trackParticle(j));
1328
1329 std::vector<const xAOD::TrackParticle*> tracksJpsi{tracksJX[0], tracksJX[1]};
1330 std::vector<const xAOD::TrackParticle*> tracksX;
1331 if(m_jxDaug_num>=3) tracksX.push_back(tracksJX[2]);
1332 if(m_jxDaug_num==4) tracksX.push_back(tracksJX[3]);
1333
1334 std::vector<double> massesV0;
1335 if(V0==LAMBDA) {
1336 massesV0 = m_massesV0_ppi;
1337 }
1338 else if(V0==LAMBDABAR) {
1339 massesV0 = m_massesV0_pip;
1340 }
1341 else if(V0==KS) {
1342 massesV0 = m_massesV0_pipi;
1343 }
1344
1345 TLorentzVector p4_moth, p4_v0, tmp;
1346 for(size_t it=0; it<JXvtx->nTrackParticles(); it++) {
1347 tmp.SetPtEtaPhiM(JXvtx->trackParticle(it)->pt(), JXvtx->trackParticle(it)->eta(), JXvtx->trackParticle(it)->phi(), massesJX[it]);
1348 p4_moth += tmp;
1349 }
1350 xAOD::BPhysHelper V0_helper(V0vtx);
1351 for(int it=0; it<V0_helper.nRefTrks(); it++) {
1352 p4_moth += V0_helper.refTrk(it,massesV0[it]);
1353 p4_v0 += V0_helper.refTrk(it,massesV0[it]);
1354 }
1355
1356 xAOD::BPhysHelper::pv_type pvtype = xAOD::BPhysHelper::PV_MIN_A0;
1357 if(m_cascadeFitWithPV==1) pvtype = xAOD::BPhysHelper::PV_MAX_SUM_PT2;
1358 else if(m_cascadeFitWithPV==2) pvtype = xAOD::BPhysHelper::PV_MIN_A0;
1359 else if(m_cascadeFitWithPV==3) pvtype = xAOD::BPhysHelper::PV_MIN_Z0;
1360 else if(m_cascadeFitWithPV==4) pvtype = xAOD::BPhysHelper::PV_MIN_Z0_BA;
1361 xAOD::BPhysHelper JX_helper(JXvtx);
1362 const xAOD::Vertex* origPv_xAOD = (m_cascadeFitWithPV>=1 && m_cascadeFitWithPV<=4) ? JX_helper.origPv(pvtype) : nullptr;
1363 const xAOD::Vertex* pv_xAOD = (m_cascadeFitWithPV>=1 && m_cascadeFitWithPV<=4) ? JX_helper.pv(pvtype) : nullptr;
1364 std::unique_ptr<Trk::RecVertex> pv_AOD;
1365 if(pv_xAOD) pv_AOD = std::make_unique<Trk::RecVertex>(pv_xAOD->position(),pv_xAOD->covariancePosition(),pv_xAOD->numberDoF(),pv_xAOD->chiSquared());
1366
1367 SG::AuxElement::Decorator<float> chi2_V1_decor("ChiSquared_V1");
1368 SG::AuxElement::Decorator<int> ndof_V1_decor("nDoF_V1");
1369 SG::AuxElement::Decorator<std::string> type_V1_decor("Type_V1");
1370
1371 SG::AuxElement::Accessor<int> mAcc_gfit("gamma_fit");
1372 SG::AuxElement::Accessor<float> mAcc_gmass("gamma_mass");
1373 SG::AuxElement::Accessor<float> mAcc_gmasserr("gamma_massError");
1374 SG::AuxElement::Accessor<float> mAcc_gchisq("gamma_chisq");
1375 SG::AuxElement::Accessor<int> mAcc_gndof("gamma_ndof");
1376 SG::AuxElement::Accessor<float> mAcc_gprob("gamma_probability");
1377
1378 SG::AuxElement::Decorator<int> mDec_gfit("gamma_fit");
1379 SG::AuxElement::Decorator<float> mDec_gmass("gamma_mass");
1380 SG::AuxElement::Decorator<float> mDec_gmasserr("gamma_massError");
1381 SG::AuxElement::Decorator<float> mDec_gchisq("gamma_chisq");
1382 SG::AuxElement::Decorator<int> mDec_gndof("gamma_ndof");
1383 SG::AuxElement::Decorator<float> mDec_gprob("gamma_probability");
1384 SG::AuxElement::Decorator< std::vector<float> > trk_pxDeco("TrackPx_V0nc");
1385 SG::AuxElement::Decorator< std::vector<float> > trk_pyDeco("TrackPy_V0nc");
1386 SG::AuxElement::Decorator< std::vector<float> > trk_pzDeco("TrackPz_V0nc");
1387
1388 std::vector<float> trk_px;
1389 std::vector<float> trk_py;
1390 std::vector<float> trk_pz;
1391
1392 if(m_extraTrk1MassHypo<=0) {
1393 double main_mass = p4_moth.M();
1394 if(m_useImprovedMass) {
1395 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_v0).M() + m_massJpsi;
1396 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_v0).M() + m_massJX;
1397 if((V0==LAMBDA || V0==LAMBDABAR) && m_massLd>0) main_mass += - p4_v0.M() + m_massLd;
1398 else if(V0==KS && m_massKs>0) main_mass += - p4_v0.M() + m_massKs;
1399 }
1400 if (main_mass < m_MassLower || main_mass > m_MassUpper) return result;
1401
1402 // Apply the user's settings to the fitter
1403 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
1404 // Robustness: http://cdsweb.cern.ch/record/685551
1405 int robustness = 0;
1406 m_iVertexFitter->setRobustness(robustness, *state);
1407 // Build up the topology
1408 // Vertex list
1409 std::vector<Trk::VertexID> vrtList;
1410 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
1411 // V0 vertex
1412 Trk::VertexID vID1;
1413 if (m_constrV0) {
1414 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,V0==KS?m_massKs:m_massLd);
1415 } else {
1416 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1417 }
1418 vrtList.push_back(vID1);
1419 Trk::VertexID vID2;
1420 if(m_JXSubVtx) {
1421 // JX vertex
1422 if (m_constrJX && m_jxDaug_num>2) {
1423 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state,m_massJX);
1424 } else {
1425 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state);
1426 }
1427 vrtList.push_back(vID2);
1428 // Mother vertex including JX and V0
1429 std::vector<const xAOD::TrackParticle*> tp;
1430 std::vector<double> tp_masses;
1431 if(m_constrMainV) {
1432 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList,*state,m_massMainV);
1433 } else {
1434 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList,*state);
1435 }
1436 }
1437 else { // m_JXSubVtx=false
1438 // Mother vertex including JX and V0
1439 if(m_constrMainV) {
1440 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state,m_massMainV);
1441 } else {
1442 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
1443 }
1444 if (m_constrJX && m_jxDaug_num>2) {
1445 std::vector<Trk::VertexID> cnstV;
1446 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
1447 ATH_MSG_WARNING("addMassConstraint for JX failed");
1448 }
1449 }
1450 }
1451 if (m_constrJpsi) {
1452 std::vector<Trk::VertexID> cnstV;
1453 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
1454 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
1455 }
1456 }
1457 if (m_constrX && m_jxDaug_num==4) {
1458 std::vector<Trk::VertexID> cnstV;
1459 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksX,cnstV,*state,m_massX).isSuccess() ) {
1460 ATH_MSG_WARNING("addMassConstraint for X failed");
1461 }
1462 }
1463 // Do the work
1464 std::unique_ptr<Trk::VxCascadeInfo> fit_result = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state, pv_AOD.get(), pv_AOD.get() && m_firstDecayAtPV ? true : false) );
1465
1466 if (fit_result) {
1467 for(auto& v : fit_result->vertices()) {
1468 if(v->nTrackParticles()==0) {
1469 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
1470 v->setTrackParticleLinks(nullLinkVector);
1471 }
1472 }
1473 // reset links to original tracks
1474 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
1475
1476 // necessary to prevent memory leak
1477 fit_result->setSVOwnership(true);
1478
1479 // Chi2/DOF cut
1480 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
1481 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
1482
1483 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
1484 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
1485 size_t iMoth = cascadeVertices.size()-1;
1486 double lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[iMoth]);
1487 if(chi2CutPassed && lxy_SV1>m_lxyV0_cut) {
1488 chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
1489 ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
1490 if(V0==LAMBDA) {
1491 type_V1_decor(*cascadeVertices[0]) = "Lambda";
1492 }
1493 else if(V0==LAMBDABAR) {
1494 type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
1495 }
1496 else if(V0==KS) {
1497 type_V1_decor(*cascadeVertices[0]) = "Ks";
1498 }
1499 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
1500 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
1501 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
1502 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
1503 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
1504 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
1505 trk_px.clear(); trk_py.clear(); trk_pz.clear();
1506 trk_px.reserve(V0_helper.nRefTrks());
1507 trk_py.reserve(V0_helper.nRefTrks());
1508 trk_pz.reserve(V0_helper.nRefTrks());
1509 for(auto&& vec3 : V0_helper.refTrks()) {
1510 trk_px.push_back( vec3.Px() );
1511 trk_py.push_back( vec3.Py() );
1512 trk_pz.push_back( vec3.Pz() );
1513 }
1514 trk_pxDeco(*cascadeVertices[0]) = trk_px;
1515 trk_pyDeco(*cascadeVertices[0]) = trk_py;
1516 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
1517
1518 result.push_back( std::make_pair(fit_result.release(),nullptr) );
1519 }
1520 }
1521 } // for m_extraTrk1MassHypo<=0
1522 else if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo<=0) {
1523 std::vector<double> massesExtra{m_extraTrk1MassHypo};
1524 std::vector<double> massesJXExtra = massesJX; massesJXExtra.push_back(m_extraTrk1MassHypo);
1525
1526 for(const xAOD::TrackParticle* tpExtra : *trackContainer) {
1527 if( tpExtra->pt()<m_extraTrk1MinPt ) continue;
1528 if( !m_trkSelector->decision(*tpExtra, nullptr) ) continue;
1529 // Check identical tracks in input
1530 if(std::find(tracksJX.cbegin(),tracksJX.cend(),tpExtra) != tracksJX.cend()) continue;
1531 if(std::find(tracksV0.cbegin(),tracksV0.cend(),tpExtra) != tracksV0.cend()) continue;
1532
1533 TLorentzVector tmp;
1534 tmp.SetPtEtaPhiM(tpExtra->pt(),tpExtra->eta(),tpExtra->phi(),m_extraTrk1MassHypo);
1535 double main_mass = (p4_moth+tmp).M();
1536 if(m_useImprovedMass) {
1537 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_v0).M() + m_massJpsi;
1538 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_v0).M() + m_massJX;
1539 if((V0==LAMBDA || V0==LAMBDABAR) && m_massLd>0) main_mass += - p4_v0.M() + m_massLd;
1540 else if(V0==KS && m_massKs>0) main_mass += - p4_v0.M() + m_massKs;
1541 }
1542 if(main_mass < m_MassLower || main_mass > m_MassUpper) continue;
1543
1544 std::vector<const xAOD::TrackParticle*> tracksExtra{tpExtra};
1545 std::vector<const xAOD::TrackParticle*> tracksJXExtra = tracksJX; tracksJXExtra.push_back(tpExtra);
1546
1547 // Apply the user's settings to the fitter
1548 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
1549 // Robustness: http://cdsweb.cern.ch/record/685551
1550 int robustness = 0;
1551 m_iVertexFitter->setRobustness(robustness, *state);
1552 // Build up the topology
1553 // Vertex list
1554 std::vector<Trk::VertexID> vrtList;
1555 std::vector<Trk::VertexID> vrtList2;
1556 Trk::VertexID vID2;
1557 if(m_JXV0SubVtx) {
1558 // V0 vertex
1559 Trk::VertexID vID1;
1560 if (m_constrV0) {
1561 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,V0==KS?m_massKs:m_massLd);
1562 } else {
1563 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1564 }
1565 vrtList.push_back(vID1);
1566 // JX+V0 vertex
1567 if(m_constrJXV0) {
1568 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state,m_massJXV0);
1569 } else {
1570 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
1571 }
1572 vrtList2.push_back(vID2);
1573 // Main vertex
1574 if(m_constrMainV) {
1575 m_iVertexFitter->nextVertex(tracksExtra,massesExtra,vrtList2,*state,m_massMainV);
1576 } else {
1577 m_iVertexFitter->nextVertex(tracksExtra,massesExtra,vrtList2,*state);
1578 }
1579 }
1580 else { // m_JXV0SubVtx==false
1581 // V0 vertex
1582 Trk::VertexID vID1;
1583 if (m_constrV0) {
1584 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,V0==KS?m_massKs:m_massLd);
1585 } else {
1586 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1587 }
1588 vrtList.push_back(vID1);
1589 if(m_JXSubVtx) {
1590 // JX vertex
1591 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state);
1592 vrtList.push_back(vID2);
1593 // Mother vertex
1594 if(m_constrMainV) {
1595 m_iVertexFitter->nextVertex(tracksExtra,massesExtra,vrtList,*state,m_massMainV);
1596 } else {
1597 m_iVertexFitter->nextVertex(tracksExtra,massesExtra,vrtList,*state);
1598 }
1599 }
1600 else { // m_JXSubVtx=false
1601 // Mother vertex includes one subvertex (V0) and JX tracks + extra track
1602 if(m_constrMainV) {
1603 vID2 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList,*state,m_massMainV);
1604 } else {
1605 vID2 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList,*state);
1606 }
1607 }
1608 }
1609 if (m_constrJX && m_jxDaug_num>2) {
1610 std::vector<Trk::VertexID> cnstV;
1611 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
1612 ATH_MSG_WARNING("addMassConstraint for JX failed");
1613 }
1614 }
1615 if (m_constrJpsi) {
1616 std::vector<Trk::VertexID> cnstV;
1617 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
1618 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
1619 }
1620 }
1621 if (m_constrX && m_jxDaug_num==4) {
1622 std::vector<Trk::VertexID> cnstV;
1623 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksX,cnstV,*state,m_massX).isSuccess() ) {
1624 ATH_MSG_WARNING("addMassConstraint for X failed");
1625 }
1626 }
1627 // Do the work
1628 std::unique_ptr<Trk::VxCascadeInfo> fit_result = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state, pv_AOD.get(), pv_AOD.get() && m_firstDecayAtPV ? true : false) );
1629
1630 if (fit_result) {
1631 for(auto& v : fit_result->vertices()) {
1632 if(v->nTrackParticles()==0) {
1633 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
1634 v->setTrackParticleLinks(nullLinkVector);
1635 }
1636 }
1637 // reset links to original tracks
1638 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
1639
1640 // necessary to prevent memory leak
1641 fit_result->setSVOwnership(true);
1642
1643 // Chi2/DOF cut
1644 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
1645 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
1646 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
1647 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
1648 size_t iMoth = cascadeVertices.size()-1;
1649 double lxy_SV1(0);
1650 if(m_JXV0SubVtx) {
1651 lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
1652 }
1653 else {
1654 lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[iMoth]);
1655 }
1656 if(chi2CutPassed && lxy_SV1>m_lxyV0_cut) {
1657 chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
1658 ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
1659 if(V0==LAMBDA) {
1660 type_V1_decor(*cascadeVertices[0]) = "Lambda";
1661 }
1662 else if(V0==LAMBDABAR) {
1663 type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
1664 }
1665 else if(V0==KS) {
1666 type_V1_decor(*cascadeVertices[0]) = "Ks";
1667 }
1668 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
1669 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
1670 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
1671 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
1672 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
1673 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
1674 trk_px.clear(); trk_py.clear(); trk_pz.clear();
1675 trk_px.reserve(V0_helper.nRefTrks());
1676 trk_py.reserve(V0_helper.nRefTrks());
1677 trk_pz.reserve(V0_helper.nRefTrks());
1678 for(auto&& vec3 : V0_helper.refTrks()) {
1679 trk_px.push_back( vec3.Px() );
1680 trk_py.push_back( vec3.Py() );
1681 trk_pz.push_back( vec3.Pz() );
1682 }
1683 trk_pxDeco(*cascadeVertices[0]) = trk_px;
1684 trk_pyDeco(*cascadeVertices[0]) = trk_py;
1685 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
1686
1687 result.push_back( std::make_pair(fit_result.release(),nullptr) );
1688 }
1689 }
1690 } // loop over trackContainer
1691 } // for m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo<=0
1692 else if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0 && m_extraTrk3MassHypo<=0) {
1693 std::vector<const xAOD::TrackParticle*> tracksPlus;
1694 std::vector<const xAOD::TrackParticle*> tracksMinus;
1695 for(const xAOD::TrackParticle* tpExtra : *trackContainer) {
1696 if( tpExtra->pt() < std::fmin(m_extraTrk1MinPt,m_extraTrk2MinPt) ) continue;
1697 if( !m_trkSelector->decision(*tpExtra, nullptr) ) continue;
1698 // Check identical tracks in input
1699 if(std::find(tracksJX.cbegin(),tracksJX.cend(),tpExtra) != tracksJX.cend()) continue;
1700 if(std::find(tracksV0.cbegin(),tracksV0.cend(),tpExtra) != tracksV0.cend()) continue;
1701 if(tpExtra->charge()>0) {
1702 tracksPlus.push_back(tpExtra);
1703 }
1704 else {
1705 tracksMinus.push_back(tpExtra);
1706 }
1707 }
1708
1709 MesonCandidateVector D0Candidates(m_maxMesonCandidates, m_MesonPtOrdering);
1710 TLorentzVector p4_ExtraTrk1, p4_ExtraTrk2;
1711 for(const xAOD::TrackParticle* tp1 : tracksPlus) {
1712 for(const xAOD::TrackParticle* tp2 : tracksMinus) {
1713 if((tp1->pt()>m_extraTrk1MinPt && tp2->pt()>m_extraTrk2MinPt) ||
1714 (tp1->pt()>m_extraTrk2MinPt && tp2->pt()>m_extraTrk1MinPt)) {
1715 p4_ExtraTrk1.SetPtEtaPhiM(tp1->pt(), tp1->eta(), tp1->phi(), m_extraTrk1MassHypo);
1716 p4_ExtraTrk2.SetPtEtaPhiM(tp2->pt(), tp2->eta(), tp2->phi(), m_extraTrk2MassHypo);
1717 double main_mass = (p4_moth+p4_ExtraTrk1+p4_ExtraTrk2).M();
1718 if(m_useImprovedMass) {
1719 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_v0).M() + m_massJpsi;
1720 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_v0).M() + m_massJX;
1721 if((V0==LAMBDA || V0==LAMBDABAR) && m_massLd>0) main_mass += - p4_v0.M() + m_massLd;
1722 else if(V0==KS && m_massKs>0) main_mass += - p4_v0.M() + m_massKs;
1723 if(m_massD0>0) main_mass += - (p4_ExtraTrk1+p4_ExtraTrk2).M() + m_massD0;
1724 }
1725 if(main_mass < m_MassLower || main_mass > m_MassUpper) continue;
1726 auto D0 = getD0Candidate(ctx,JXvtx,tp1,tp2);
1727 if(D0.extraTrack1) D0Candidates.push_back(D0);
1728 }
1729 }
1730 }
1731
1732 std::vector<double> massesExtra{m_extraTrk1MassHypo,m_extraTrk2MassHypo};
1733
1734 for(auto&& D0 : D0Candidates.vector()) {
1735 std::vector<const xAOD::TrackParticle*> tracksExtra{D0.extraTrack1,D0.extraTrack2};
1736
1737 // Apply the user's settings to the fitter
1738 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
1739 // Robustness: http://cdsweb.cern.ch/record/685551
1740 int robustness = 0;
1741 m_iVertexFitter->setRobustness(robustness, *state);
1742 // Build up the topology
1743 // Vertex list
1744 std::vector<Trk::VertexID> vrtList;
1745 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
1746 // V0 vertex
1747 Trk::VertexID vID1;
1748 if (m_constrV0) {
1749 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,V0==KS?m_massKs:m_massLd);
1750 } else {
1751 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1752 }
1753 vrtList.push_back(vID1);
1754 // D0 vertex
1755 Trk::VertexID vID2;
1756 if (m_constrD0) {
1757 vID2 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state,m_massD0);
1758 } else {
1759 vID2 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state);
1760 }
1761 vrtList.push_back(vID2);
1762 // Mother vertex includes one subvertex (V0) and JX tracks + extra track
1763 Trk::VertexID vID3;
1764 if(m_constrMainV) {
1765 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state,m_massMainV);
1766 } else {
1767 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
1768 }
1769 if (m_constrJX && m_jxDaug_num>2) {
1770 std::vector<Trk::VertexID> cnstV;
1771 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
1772 ATH_MSG_WARNING("addMassConstraint for JX failed");
1773 }
1774 }
1775 if (m_constrJpsi) {
1776 std::vector<Trk::VertexID> cnstV;
1777 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
1778 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
1779 }
1780 }
1781 if (m_constrX && m_jxDaug_num==4) {
1782 std::vector<Trk::VertexID> cnstV;
1783 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksX,cnstV,*state,m_massX).isSuccess() ) {
1784 ATH_MSG_WARNING("addMassConstraint for X failed");
1785 }
1786 }
1787 // Do the work
1788 std::unique_ptr<Trk::VxCascadeInfo> fit_result = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state, pv_AOD.get(), pv_AOD.get() && m_firstDecayAtPV ? true : false) );
1789
1790 if (fit_result) {
1791 for(auto& v : fit_result->vertices()) {
1792 if(v->nTrackParticles()==0) {
1793 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
1794 v->setTrackParticleLinks(nullLinkVector);
1795 }
1796 }
1797 // reset links to original tracks
1798 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
1799
1800 // necessary to prevent memory leak
1801 fit_result->setSVOwnership(true);
1802
1803 // Chi2/DOF cut
1804 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
1805 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
1806 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
1807 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
1808 size_t iMoth = cascadeVertices.size()-1;
1809 double lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[iMoth]);
1810 double lxy_SV2 = m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[iMoth]);
1811 if(chi2CutPassed && lxy_SV1>m_lxyV0_cut && lxy_SV2>m_lxyD0_cut) {
1812 chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
1813 ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
1814 if(V0==LAMBDA) {
1815 type_V1_decor(*cascadeVertices[0]) = "Lambda";
1816 }
1817 else if(V0==LAMBDABAR) {
1818 type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
1819 }
1820 else if(V0==KS) {
1821 type_V1_decor(*cascadeVertices[0]) = "Ks";
1822 }
1823 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
1824 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
1825 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
1826 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
1827 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
1828 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
1829 trk_px.clear(); trk_py.clear(); trk_pz.clear();
1830 trk_px.reserve(V0_helper.nRefTrks());
1831 trk_py.reserve(V0_helper.nRefTrks());
1832 trk_pz.reserve(V0_helper.nRefTrks());
1833 for(auto&& vec3 : V0_helper.refTrks()) {
1834 trk_px.push_back( vec3.Px() );
1835 trk_py.push_back( vec3.Py() );
1836 trk_pz.push_back( vec3.Pz() );
1837 }
1838 trk_pxDeco(*cascadeVertices[0]) = trk_px;
1839 trk_pyDeco(*cascadeVertices[0]) = trk_py;
1840 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
1841
1842 result.push_back( std::make_pair(fit_result.release(),nullptr) );
1843 }
1844 }
1845 } // loop over D0Candidates
1846 } // for m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0 && m_extraTrk3MassHypo<=0
1847 else if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0 && m_extraTrk3MassHypo>0) {
1848 std::vector<const xAOD::TrackParticle*> tracksPlus;
1849 std::vector<const xAOD::TrackParticle*> tracksMinus;
1850 double minTrkPt = std::fmin(std::fmin(m_extraTrk1MinPt,m_extraTrk2MinPt),m_extraTrk3MinPt);
1851 for(const xAOD::TrackParticle* tpExtra : *trackContainer) {
1852 if( tpExtra->pt() < minTrkPt ) continue;
1853 if( !m_trkSelector->decision(*tpExtra, nullptr) ) continue;
1854 // Check identical tracks in input
1855 if(std::find(tracksJX.cbegin(),tracksJX.cend(),tpExtra) != tracksJX.cend()) continue;
1856 if(std::find(tracksV0.cbegin(),tracksV0.cend(),tpExtra) != tracksV0.cend()) continue;
1857 if(tpExtra->charge()>0) {
1858 tracksPlus.push_back(tpExtra);
1859 }
1860 else {
1861 tracksMinus.push_back(tpExtra);
1862 }
1863 }
1864
1865 MesonCandidateVector DpmCandidates(m_maxMesonCandidates, m_MesonPtOrdering);
1866 TLorentzVector p4_ExtraTrk1, p4_ExtraTrk2, p4_ExtraTrk3;
1867 // +,-,- combination (D- -> K+ pi- pi-)
1868 for(const xAOD::TrackParticle* tp1 : tracksPlus) {
1869 if( tp1->pt() < m_extraTrk1MinPt ) continue;
1870 for(auto tp2Itr=tracksMinus.cbegin(); tp2Itr!=tracksMinus.cend(); ++tp2Itr) {
1871 const xAOD::TrackParticle* tp2 = *tp2Itr;
1872 for(auto tp3Itr=tp2Itr+1; tp3Itr!=tracksMinus.cend(); ++tp3Itr) {
1873 const xAOD::TrackParticle* tp3 = *tp3Itr;
1874 if((tp2->pt()>m_extraTrk2MinPt && tp3->pt()>m_extraTrk3MinPt) ||
1875 (tp2->pt()>m_extraTrk3MinPt && tp3->pt()>m_extraTrk2MinPt)) {
1876 p4_ExtraTrk1.SetPtEtaPhiM(tp1->pt(), tp1->eta(), tp1->phi(), m_extraTrk1MassHypo);
1877 p4_ExtraTrk2.SetPtEtaPhiM(tp2->pt(), tp2->eta(), tp2->phi(), m_extraTrk2MassHypo);
1878 p4_ExtraTrk3.SetPtEtaPhiM(tp3->pt(), tp3->eta(), tp3->phi(), m_extraTrk3MassHypo);
1879 double main_mass = (p4_moth+p4_ExtraTrk1+p4_ExtraTrk2+p4_ExtraTrk3).M();
1880 if(m_useImprovedMass) {
1881 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_v0).M() + m_massJpsi;
1882 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_v0).M() + m_massJX;
1883 if((V0==LAMBDA || V0==LAMBDABAR) && m_massLd>0) main_mass += - p4_v0.M() + m_massLd;
1884 else if(V0==KS && m_massKs>0) main_mass += - p4_v0.M() + m_massKs;
1885 if(m_massDpm>0) main_mass += - (p4_ExtraTrk1+p4_ExtraTrk2+p4_ExtraTrk3).M() + m_massDpm;
1886 }
1887 if(main_mass < m_MassLower || main_mass > m_MassUpper) continue;
1888 auto Dpm = getDpmCandidate(ctx,JXvtx,tp1,tp2,tp3);
1889 if(Dpm.extraTrack1) DpmCandidates.push_back(Dpm);
1890 }
1891 }
1892 }
1893 }
1894 // -,+,+ combination (D+ -> K- pi+ pi+)
1895 for(const xAOD::TrackParticle* tp1 : tracksMinus) {
1896 if( tp1->pt() < m_extraTrk1MinPt ) continue;
1897 for(auto tp2Itr=tracksPlus.cbegin(); tp2Itr!=tracksPlus.cend(); ++tp2Itr) {
1898 const xAOD::TrackParticle* tp2 = *tp2Itr;
1899 for(auto tp3Itr=tp2Itr+1; tp3Itr!=tracksPlus.cend(); ++tp3Itr) {
1900 const xAOD::TrackParticle* tp3 = *tp3Itr;
1901 if((tp2->pt()>m_extraTrk2MinPt && tp3->pt()>m_extraTrk3MinPt) ||
1902 (tp2->pt()>m_extraTrk3MinPt && tp3->pt()>m_extraTrk2MinPt)) {
1903 p4_ExtraTrk1.SetPtEtaPhiM(tp1->pt(), tp1->eta(), tp1->phi(), m_extraTrk1MassHypo);
1904 p4_ExtraTrk2.SetPtEtaPhiM(tp2->pt(), tp2->eta(), tp2->phi(), m_extraTrk2MassHypo);
1905 p4_ExtraTrk3.SetPtEtaPhiM(tp3->pt(), tp3->eta(), tp3->phi(), m_extraTrk3MassHypo);
1906 double main_mass = (p4_moth+p4_ExtraTrk1+p4_ExtraTrk2+p4_ExtraTrk3).M();
1907 if(m_useImprovedMass) {
1908 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_v0).M() + m_massJpsi;
1909 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_v0).M() + m_massJX;
1910 if((V0==LAMBDA || V0==LAMBDABAR) && m_massLd>0) main_mass += - p4_v0.M() + m_massLd;
1911 else if(V0==KS && m_massKs>0) main_mass += - p4_v0.M() + m_massKs;
1912 if(m_massDpm>0) main_mass += - (p4_ExtraTrk1+p4_ExtraTrk2+p4_ExtraTrk3).M() + m_massDpm;
1913 }
1914 if(main_mass < m_MassLower || main_mass > m_MassUpper) continue;
1915 auto Dpm = getDpmCandidate(ctx,JXvtx,tp1,tp2,tp3);
1916 if(Dpm.extraTrack1) DpmCandidates.push_back(Dpm);
1917 }
1918 }
1919 }
1920 }
1921
1922 std::vector<double> massesExtra{m_extraTrk1MassHypo,m_extraTrk2MassHypo,m_extraTrk3MassHypo};
1923
1924 for(auto&& Dpm : DpmCandidates.vector()) {
1925 std::vector<const xAOD::TrackParticle*> tracksExtra{Dpm.extraTrack1,Dpm.extraTrack2,Dpm.extraTrack3};
1926
1927 // Apply the user's settings to the fitter
1928 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
1929 // Robustness: http://cdsweb.cern.ch/record/685551
1930 int robustness = 0;
1931 m_iVertexFitter->setRobustness(robustness, *state);
1932 // Build up the topology
1933 // Vertex list
1934 std::vector<Trk::VertexID> vrtList;
1935 std::vector<Trk::VertexID> vrtList2;
1936 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
1937 if(m_JXV0SubVtx) {
1938 // V0 vertex
1939 Trk::VertexID vID1;
1940 if (m_constrV0) {
1941 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,V0==KS?m_massKs:m_massLd);
1942 } else {
1943 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1944 }
1945 vrtList.push_back(vID1);
1946 // JX+V0 vertex
1947 Trk::VertexID vID2;
1948 if(m_constrJXV0) {
1949 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state,m_massJXV0);
1950 } else {
1951 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
1952 }
1953 vrtList2.push_back(vID2);
1954 Trk::VertexID vID3;
1955 if (m_constrDpm) {
1956 vID3 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state,m_massDpm);
1957 } else {
1958 vID3 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state);
1959 }
1960 vrtList2.push_back(vID3);
1961 // Mother vertex includes two subvertices: V0, Dpm and JX tracks
1962 std::vector<const xAOD::TrackParticle*> tp;
1963 std::vector<double> tp_masses;
1964 if(m_constrMainV) {
1965 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2,*state,m_massMainV);
1966 } else {
1967 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2,*state);
1968 }
1969 if (m_constrJX && m_jxDaug_num>2) {
1970 std::vector<Trk::VertexID> cnstV;
1971 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
1972 ATH_MSG_WARNING("addMassConstraint for JX failed");
1973 }
1974 }
1975 if (m_constrJpsi) {
1976 std::vector<Trk::VertexID> cnstV;
1977 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
1978 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
1979 }
1980 }
1981 if (m_constrX && m_jxDaug_num==4) {
1982 std::vector<Trk::VertexID> cnstV;
1983 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksX,cnstV,*state,m_massX).isSuccess() ) {
1984 ATH_MSG_WARNING("addMassConstraint for X failed");
1985 }
1986 }
1987 }
1988 else { // m_JXV0SubVtx==false
1989 // V0 vertex
1990 Trk::VertexID vID1;
1991 if (m_constrV0) {
1992 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,V0==KS?m_massKs:m_massLd);
1993 } else {
1994 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1995 }
1996 vrtList.push_back(vID1);
1997 // Dpm vertex
1998 Trk::VertexID vID2;
1999 if (m_constrDpm) {
2000 vID2 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state,m_massDpm);
2001 } else {
2002 vID2 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state);
2003 }
2004 vrtList.push_back(vID2);
2005 // Mother vertex includes two subvertices: V0, Dpm and JX tracks
2006 Trk::VertexID vID3;
2007 if(m_constrMainV) {
2008 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state,m_massMainV);
2009 } else {
2010 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
2011 }
2012 if (m_constrJX && m_jxDaug_num>2) {
2013 std::vector<Trk::VertexID> cnstV;
2014 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
2015 ATH_MSG_WARNING("addMassConstraint for JX failed");
2016 }
2017 }
2018 if (m_constrJpsi) {
2019 std::vector<Trk::VertexID> cnstV;
2020 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
2021 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
2022 }
2023 }
2024 if (m_constrX && m_jxDaug_num==4) {
2025 std::vector<Trk::VertexID> cnstV;
2026 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksX,cnstV,*state,m_massX).isSuccess() ) {
2027 ATH_MSG_WARNING("addMassConstraint for X failed");
2028 }
2029 }
2030 }
2031 // Do the work
2032 std::unique_ptr<Trk::VxCascadeInfo> fit_result = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state, pv_AOD.get(), pv_AOD.get() && m_firstDecayAtPV ? true : false) );
2033
2034 if (fit_result) {
2035 for(auto& v : fit_result->vertices()) {
2036 if(v->nTrackParticles()==0) {
2037 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
2038 v->setTrackParticleLinks(nullLinkVector);
2039 }
2040 }
2041 // reset links to original tracks
2042 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
2043
2044 // necessary to prevent memory leak
2045 fit_result->setSVOwnership(true);
2046
2047 // Chi2/DOF cut
2048 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
2049 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
2050 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
2051 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
2052 size_t iMoth = cascadeVertices.size()-1;
2053 double lxy_SV1(0), lxy_SV2(0);
2054 if(m_JXV0SubVtx) {
2055 lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
2056 lxy_SV2 = m_CascadeTools->lxy(moms[2],cascadeVertices[2],cascadeVertices[iMoth]);
2057 }
2058 else {
2059 lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[iMoth]);
2060 lxy_SV2 = m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[iMoth]);
2061 }
2062 if(chi2CutPassed && lxy_SV1>m_lxyV0_cut && lxy_SV2>m_lxyDpm_cut) {
2063 chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
2064 ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
2065 if(V0==LAMBDA) {
2066 type_V1_decor(*cascadeVertices[0]) = "Lambda";
2067 }
2068 else if(V0==LAMBDABAR) {
2069 type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
2070 }
2071 else if(V0==KS) {
2072 type_V1_decor(*cascadeVertices[0]) = "Ks";
2073 }
2074 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
2075 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
2076 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
2077 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
2078 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
2079 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
2080 trk_px.clear(); trk_py.clear(); trk_pz.clear();
2081 trk_px.reserve(V0_helper.nRefTrks());
2082 trk_py.reserve(V0_helper.nRefTrks());
2083 trk_pz.reserve(V0_helper.nRefTrks());
2084 for(auto&& vec3 : V0_helper.refTrks()) {
2085 trk_px.push_back( vec3.Px() );
2086 trk_py.push_back( vec3.Py() );
2087 trk_pz.push_back( vec3.Pz() );
2088 }
2089 trk_pxDeco(*cascadeVertices[0]) = trk_px;
2090 trk_pyDeco(*cascadeVertices[0]) = trk_py;
2091 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
2092
2093 // refit with main vertex mass constraint if required
2094 if(!m_constrMainV && m_doPostMainVContrFit) {
2095 TLorentzVector totalMom;
2096 for(size_t it=0; it<moms[iMoth].size(); it++) totalMom += moms[iMoth][it];
2097 double mainV_mass = totalMom.M();
2098 if(mainV_mass>m_PostMassLower && mainV_mass<m_PostMassUpper) {
2099 std::unique_ptr<Trk::IVKalState> state_mvc = m_iVertexFitter->makeState(ctx);
2100 int robustness_mvc = 0;
2101 m_iVertexFitter->setRobustness(robustness_mvc, *state_mvc);
2102 std::vector<Trk::VertexID> vrtList_mvc;
2103 std::vector<Trk::VertexID> vrtList2_mvc;
2104 if(m_JXV0SubVtx) {
2105 Trk::VertexID vID1_mvc; // V0 vertex
2106 if (m_constrV0) {
2107 vID1_mvc = m_iVertexFitter->startVertex(tracksV0,massesV0,*state_mvc,V0==KS?m_massKs:m_massLd);
2108 } else {
2109 vID1_mvc = m_iVertexFitter->startVertex(tracksV0,massesV0,*state_mvc);
2110 }
2111 vrtList_mvc.push_back(vID1_mvc);
2112 Trk::VertexID vID2_mvc; // JX+V0 vertex
2113 if(m_constrJXV0) {
2114 vID2_mvc = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList_mvc,*state_mvc,m_massJXV0);
2115 } else {
2116 vID2_mvc = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList_mvc,*state_mvc);
2117 }
2118 vrtList2_mvc.push_back(vID2_mvc);
2119 Trk::VertexID vID3_mvc; // Dpm vertex
2120 if (m_constrDpm) {
2121 vID3_mvc = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state_mvc,m_massDpm);
2122 } else {
2123 vID3_mvc = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state_mvc);
2124 }
2125 vrtList2_mvc.push_back(vID3_mvc);
2126 std::vector<const xAOD::TrackParticle*> tp;
2127 std::vector<double> tp_masses;
2128 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2_mvc,*state_mvc,m_massMainV); // Mother vertex
2129 if (m_constrJX && m_jxDaug_num>2) {
2130 std::vector<Trk::VertexID> cnstV_mvc;
2131 if ( !m_iVertexFitter->addMassConstraint(vID2_mvc,tracksJX,cnstV_mvc,*state_mvc,m_massJX).isSuccess() ) {
2132 ATH_MSG_WARNING("addMassConstraint for JX failed");
2133 }
2134 }
2135 if (m_constrJpsi) {
2136 std::vector<Trk::VertexID> cnstV_mvc;
2137 if ( !m_iVertexFitter->addMassConstraint(vID2_mvc,tracksJpsi,cnstV_mvc,*state_mvc,m_massJpsi).isSuccess() ) {
2138 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
2139 }
2140 }
2141 if (m_constrX && m_jxDaug_num==4) {
2142 std::vector<Trk::VertexID> cnstV_mvc;
2143 if ( !m_iVertexFitter->addMassConstraint(vID2_mvc,tracksX,cnstV_mvc,*state_mvc,m_massX).isSuccess() ) {
2144 ATH_MSG_WARNING("addMassConstraint for X failed");
2145 }
2146 }
2147 }
2148 else {
2149 Trk::VertexID vID1_mvc; // V0 vertex
2150 if (m_constrV0) {
2151 vID1_mvc = m_iVertexFitter->startVertex(tracksV0,massesV0,*state_mvc,V0==KS?m_massKs:m_massLd);
2152 } else {
2153 vID1_mvc = m_iVertexFitter->startVertex(tracksV0,massesV0,*state_mvc);
2154 }
2155 vrtList_mvc.push_back(vID1_mvc);
2156 Trk::VertexID vID2_mvc; // Dpm vertex
2157 if (m_constrDpm) {
2158 vID2_mvc = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state_mvc,m_massDpm);
2159 } else {
2160 vID2_mvc = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state_mvc);
2161 }
2162 vrtList_mvc.push_back(vID2_mvc);
2163 Trk::VertexID vID3_mvc = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList_mvc,*state_mvc,m_massMainV); // Mother vertex
2164 if (m_constrJX && m_jxDaug_num>2) {
2165 std::vector<Trk::VertexID> cnstV_mvc;
2166 if ( !m_iVertexFitter->addMassConstraint(vID3_mvc,tracksJX,cnstV_mvc,*state_mvc,m_massJX).isSuccess() ) {
2167 ATH_MSG_WARNING("addMassConstraint for JX failed");
2168 }
2169 }
2170 if (m_constrJpsi) {
2171 std::vector<Trk::VertexID> cnstV_mvc;
2172 if ( !m_iVertexFitter->addMassConstraint(vID3_mvc,tracksJpsi,cnstV_mvc,*state_mvc,m_massJpsi).isSuccess() ) {
2173 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
2174 }
2175 }
2176 if (m_constrX && m_jxDaug_num==4) {
2177 std::vector<Trk::VertexID> cnstV_mvc;
2178 if ( !m_iVertexFitter->addMassConstraint(vID3_mvc,tracksX,cnstV_mvc,*state_mvc,m_massX).isSuccess() ) {
2179 ATH_MSG_WARNING("addMassConstraint for X failed");
2180 }
2181 }
2182 }
2183 // Do the work
2184 std::unique_ptr<Trk::VxCascadeInfo> fit_result_mvc = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state_mvc, pv_AOD.get(), pv_AOD.get() && m_firstDecayAtPV ? true : false) );
2185
2186 if (fit_result_mvc) {
2187 for(auto& v : fit_result_mvc->vertices()) {
2188 if(v->nTrackParticles()==0) {
2189 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
2190 v->setTrackParticleLinks(nullLinkVector);
2191 }
2192 }
2193 BPhysPVCascadeTools::PrepareVertexLinks(fit_result_mvc.get(), trackCols);
2194 fit_result_mvc->setSVOwnership(true);
2195 chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
2196 ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
2197 if(V0==LAMBDA) {
2198 type_V1_decor(*cascadeVertices[0]) = "Lambda";
2199 }
2200 else if(V0==LAMBDABAR) {
2201 type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
2202 }
2203 else if(V0==KS) {
2204 type_V1_decor(*cascadeVertices[0]) = "Ks";
2205 }
2206 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
2207 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
2208 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
2209 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
2210 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
2211 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
2212 trk_px.clear(); trk_py.clear(); trk_pz.clear();
2213 trk_px.reserve(V0_helper.nRefTrks());
2214 trk_py.reserve(V0_helper.nRefTrks());
2215 trk_pz.reserve(V0_helper.nRefTrks());
2216 for(auto&& vec3 : V0_helper.refTrks()) {
2217 trk_px.push_back( vec3.Px() );
2218 trk_py.push_back( vec3.Py() );
2219 trk_pz.push_back( vec3.Pz() );
2220 }
2221 trk_pxDeco(*cascadeVertices[0]) = trk_px;
2222 trk_pyDeco(*cascadeVertices[0]) = trk_py;
2223 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
2224
2225 result.push_back( std::make_pair(fit_result.release(),fit_result_mvc.release()) );
2226 }
2227 else result.push_back( std::make_pair(fit_result.release(),nullptr) );
2228 }
2229 else result.push_back( std::make_pair(fit_result.release(),nullptr) );
2230 }
2231 else result.push_back( std::make_pair(fit_result.release(),nullptr) );
2232 }
2233 }
2234 } // loop over DpmCandidates
2235 } // for m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0 && m_extraTrk3MassHypo>0
2236
2237 if(pv_xAOD) {
2238 for(auto cascade_info_pair : result) {
2239 if(cascade_info_pair.first && cascade_info_pair.first->getParticleMoms().size()>0) {
2240 size_t index = cascade_info_pair.first->getParticleMoms().size() - 1;
2241 const std::vector<TLorentzVector> &mom = cascade_info_pair.first->getParticleMoms()[index];
2242 const Amg::MatrixX &cov = cascade_info_pair.first->getCovariance()[index];
2243 const xAOD::Vertex* mainVertex = cascade_info_pair.first->vertices()[index];
2244 xAOD::BPhysHypoHelper vtx(m_hypoName, mainVertex);
2245 bool isInDefaultPVCont = false;
2246 for(const xAOD::Vertex* pvVtx : *defaultPVContainer) {
2247 if(pv_xAOD == pvVtx) { isInDefaultPVCont = true; break; }
2248 }
2249 if(isInDefaultPVCont) vtx.setPv( pv_xAOD, defaultPVContainer, pvtype );
2250 else vtx.setPv( pv_xAOD, pvContainer, pvtype );
2251 if(origPv_xAOD) vtx.setOrigPv( origPv_xAOD, defaultPVContainer, pvtype );
2252 vtx.setLxy ( m_CascadeTools->lxy (mom, vtx.vtx(), pv_xAOD), pvtype );
2253 vtx.setLxyErr ( m_CascadeTools->lxyError (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2254 vtx.setA0 ( m_CascadeTools->a0 (mom, vtx.vtx(), pv_xAOD), pvtype );
2255 vtx.setA0Err ( m_CascadeTools->a0Error (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2256 vtx.setA0xy ( m_CascadeTools->a0xy (mom, vtx.vtx(), pv_xAOD), pvtype );
2257 vtx.setA0xyErr( m_CascadeTools->a0xyError(mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2258 vtx.setZ0 ( m_CascadeTools->a0z (mom, vtx.vtx(), pv_xAOD), pvtype );
2259 vtx.setZ0Err ( m_CascadeTools->a0zError (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2260 vtx.setRefitPVStatus( 0, pvtype );
2261 // Proper decay times
2262 vtx.setTau( m_CascadeTools->tau(mom, vtx.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
2263 vtx.setTauErr( m_CascadeTools->tauError(mom, cov, vtx.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
2264 vtx.setTau( m_CascadeTools->tau(mom, vtx.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
2265 vtx.setTauErr( m_CascadeTools->tauError(mom, cov, vtx.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
2266
2267 if(cascade_info_pair.second && cascade_info_pair.second->getParticleMoms().size()>0) {
2268 index = cascade_info_pair.second->getParticleMoms().size() - 1;
2269 const std::vector<TLorentzVector> &mom_mvc = cascade_info_pair.second->getParticleMoms()[index];
2270 const Amg::MatrixX &cov_mvc = cascade_info_pair.second->getCovariance()[index];
2271 const xAOD::Vertex* mainVertex_mvc = cascade_info_pair.second->vertices()[index];
2272 xAOD::BPhysHypoHelper vtx_mvc(m_hypoName, mainVertex_mvc);
2273 if(isInDefaultPVCont) vtx_mvc.setPv( pv_xAOD, defaultPVContainer, pvtype );
2274 else vtx_mvc.setPv( pv_xAOD, pvContainer, pvtype );
2275 if(origPv_xAOD) vtx.setOrigPv( origPv_xAOD, defaultPVContainer, pvtype );
2276 vtx_mvc.setLxy ( m_CascadeTools->lxy (mom_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2277 vtx_mvc.setLxyErr ( m_CascadeTools->lxyError (mom_mvc, cov_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2278 vtx_mvc.setA0 ( m_CascadeTools->a0 (mom_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2279 vtx_mvc.setA0Err ( m_CascadeTools->a0Error (mom_mvc, cov_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2280 vtx_mvc.setA0xy ( m_CascadeTools->a0xy (mom_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2281 vtx_mvc.setA0xyErr( m_CascadeTools->a0xyError(mom_mvc, cov_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2282 vtx_mvc.setZ0 ( m_CascadeTools->a0z (mom_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2283 vtx_mvc.setZ0Err ( m_CascadeTools->a0zError (mom_mvc, cov_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2284 vtx_mvc.setRefitPVStatus( 0, pvtype );
2285 // Proper decay times
2286 vtx_mvc.setTau( m_CascadeTools->tau(mom_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
2287 vtx_mvc.setTauErr( m_CascadeTools->tauError(mom_mvc, cov_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
2288 vtx_mvc.setTau( m_CascadeTools->tau(mom_mvc, vtx_mvc.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
2289 vtx_mvc.setTauErr( m_CascadeTools->tauError(mom_mvc, cov_mvc, vtx_mvc.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
2290 }
2291 }
2292 }
2293 }
2294
2295 return result;
2296 }
2297
2298 std::vector<std::pair<Trk::VxCascadeInfo*,Trk::VxCascadeInfo*> > JpsiXPlusDisplaced::fitMainVtx(const EventContext& ctx, const xAOD::Vertex* JXvtx, const std::vector<double>& massesJX, const XiCandidate& disVtx, const xAOD::TrackParticleContainer* trackContainer, const std::vector<const xAOD::TrackParticleContainer*>& trackCols, const xAOD::VertexContainer* defaultPVContainer, const xAOD::VertexContainer* pvContainer) const {
2299 std::vector<std::pair<Trk::VxCascadeInfo*,Trk::VxCascadeInfo*> > result;
2300
2301 std::vector<const xAOD::TrackParticle*> tracksJX;
2302 tracksJX.reserve(JXvtx->nTrackParticles());
2303 for(size_t i=0; i<JXvtx->nTrackParticles(); i++) tracksJX.push_back(JXvtx->trackParticle(i));
2304 if (tracksJX.size() != massesJX.size()) {
2305 ATH_MSG_ERROR("Problems with JX input: number of tracks or track mass inputs is not correct!");
2306 return result;
2307 }
2308 // Check identical tracks in input
2309 if(std::find(tracksJX.cbegin(), tracksJX.cend(), disVtx.V0vtx->trackParticle(0)) != tracksJX.cend()) return result;
2310 if(std::find(tracksJX.cbegin(), tracksJX.cend(), disVtx.V0vtx->trackParticle(1)) != tracksJX.cend()) return result;
2311 std::vector<const xAOD::TrackParticle*> tracksV0;
2312 tracksV0.reserve(disVtx.V0vtx->nTrackParticles());
2313 for(size_t j=0; j<disVtx.V0vtx->nTrackParticles(); j++) tracksV0.push_back(disVtx.V0vtx->trackParticle(j));
2314
2315 if(std::find(tracksJX.cbegin(), tracksJX.cend(), disVtx.track) != tracksJX.cend()) return result;
2316 std::vector<const xAOD::TrackParticle*> tracks3{disVtx.track};
2317 std::vector<double> massesDis3{m_disVDaug3MassHypo};
2318
2319 std::vector<const xAOD::TrackParticle*> tracksJpsi{tracksJX[0], tracksJX[1]};
2320 std::vector<const xAOD::TrackParticle*> tracksX;
2321 if(m_jxDaug_num>=3) tracksX.push_back(tracksJX[2]);
2322 if(m_jxDaug_num==4) tracksX.push_back(tracksJX[3]);
2323
2324 std::vector<double> massesV0;
2325 if(disVtx.V0type==LAMBDA) {
2326 massesV0 = m_massesV0_ppi;
2327 }
2328 else if(disVtx.V0type==LAMBDABAR) {
2329 massesV0 = m_massesV0_pip;
2330 }
2331 else if(disVtx.V0type==KS) {
2332 massesV0 = m_massesV0_pipi;
2333 }
2334
2335 std::vector<double> massesDisV = massesV0; massesDisV.push_back(m_disVDaug3MassHypo);
2336
2337 TLorentzVector p4_moth, p4_disV, tmp;
2338 for(size_t it=0; it<JXvtx->nTrackParticles(); it++) {
2339 tmp.SetPtEtaPhiM(JXvtx->trackParticle(it)->pt(), JXvtx->trackParticle(it)->eta(), JXvtx->trackParticle(it)->phi(), massesJX[it]);
2340 p4_moth += tmp;
2341 }
2342 p4_disV += disVtx.p4_V0track1; p4_disV += disVtx.p4_V0track2; p4_disV += disVtx.p4_disVtrack;
2343 p4_moth += p4_disV;
2344
2345 xAOD::BPhysHelper::pv_type pvtype = xAOD::BPhysHelper::PV_MIN_A0;
2346 if(m_cascadeFitWithPV==1) pvtype = xAOD::BPhysHelper::PV_MAX_SUM_PT2;
2347 else if(m_cascadeFitWithPV==2) pvtype = xAOD::BPhysHelper::PV_MIN_A0;
2348 else if(m_cascadeFitWithPV==3) pvtype = xAOD::BPhysHelper::PV_MIN_Z0;
2349 else if(m_cascadeFitWithPV==4) pvtype = xAOD::BPhysHelper::PV_MIN_Z0_BA;
2350 xAOD::BPhysHelper JX_helper(JXvtx);
2351 const xAOD::Vertex* origPv_xAOD = (m_cascadeFitWithPV>=1 && m_cascadeFitWithPV<=4) ? JX_helper.origPv(pvtype) : nullptr;
2352 const xAOD::Vertex* pv_xAOD = (m_cascadeFitWithPV>=1 && m_cascadeFitWithPV<=4) ? JX_helper.pv(pvtype) : nullptr;
2353 std::unique_ptr<Trk::RecVertex> pv_AOD;
2354 if(pv_xAOD) pv_AOD = std::make_unique<Trk::RecVertex>(pv_xAOD->position(),pv_xAOD->covariancePosition(),pv_xAOD->numberDoF(),pv_xAOD->chiSquared());
2355
2356 SG::AuxElement::Decorator<float> chi2_V0_decor("ChiSquared_V0");
2357 SG::AuxElement::Decorator<int> ndof_V0_decor("nDoF_V0");
2358 SG::AuxElement::Decorator<std::string> type_V0_decor("Type_V0");
2359
2360 SG::AuxElement::Accessor<int> mAcc_gfit("gamma_fit");
2361 SG::AuxElement::Accessor<float> mAcc_gmass("gamma_mass");
2362 SG::AuxElement::Accessor<float> mAcc_gmasserr("gamma_massError");
2363 SG::AuxElement::Accessor<float> mAcc_gchisq("gamma_chisq");
2364 SG::AuxElement::Accessor<int> mAcc_gndof("gamma_ndof");
2365 SG::AuxElement::Accessor<float> mAcc_gprob("gamma_probability");
2366
2367 SG::AuxElement::Decorator<int> mDec_gfit("gamma_fit");
2368 SG::AuxElement::Decorator<float> mDec_gmass("gamma_mass");
2369 SG::AuxElement::Decorator<float> mDec_gmasserr("gamma_massError");
2370 SG::AuxElement::Decorator<float> mDec_gchisq("gamma_chisq");
2371 SG::AuxElement::Decorator<int> mDec_gndof("gamma_ndof");
2372 SG::AuxElement::Decorator<float> mDec_gprob("gamma_probability");
2373 SG::AuxElement::Decorator< std::vector<float> > trk_pxDeco("TrackPx_V0nc");
2374 SG::AuxElement::Decorator< std::vector<float> > trk_pyDeco("TrackPy_V0nc");
2375 SG::AuxElement::Decorator< std::vector<float> > trk_pzDeco("TrackPz_V0nc");
2376 SG::AuxElement::Decorator<float> trk_px_deco("TrackPx_DisVnc");
2377 SG::AuxElement::Decorator<float> trk_py_deco("TrackPy_DisVnc");
2378 SG::AuxElement::Decorator<float> trk_pz_deco("TrackPz_DisVnc");
2379
2380 std::vector<float> trk_px;
2381 std::vector<float> trk_py;
2382 std::vector<float> trk_pz;
2383
2384 if(m_extraTrk1MassHypo<=0) {
2385 double main_mass = p4_moth.M();
2386 if(m_useImprovedMass) {
2387 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_disV).M() + m_massJpsi;
2388 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_disV).M() + m_massJX;
2389 if(m_massDisV>0) main_mass += - p4_disV.M() + m_massDisV;
2390 }
2391 if (main_mass < m_MassLower || main_mass > m_MassUpper) return result;
2392
2393 // Apply the user's settings to the fitter
2394 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
2395 // Robustness: http://cdsweb.cern.ch/record/685551
2396 int robustness = 0;
2397 m_iVertexFitter->setRobustness(robustness, *state);
2398 // Build up the topology
2399 // Vertex list
2400 std::vector<Trk::VertexID> vrtList;
2401 std::vector<Trk::VertexID> vrtList2;
2402 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
2403 // V0 vertex
2404 Trk::VertexID vID1;
2405 if (m_constrV0) {
2406 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,disVtx.V0type==KS?m_massKs:m_massLd);
2407 } else {
2408 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
2409 }
2410 vrtList.push_back(vID1);
2411 // Displaced vertex
2412 Trk::VertexID vID2;
2413 if (m_constrDisV) {
2414 vID2 = m_iVertexFitter->nextVertex(tracks3,massesDis3,vrtList,*state,m_massDisV);
2415 } else {
2416 vID2 = m_iVertexFitter->nextVertex(tracks3,massesDis3,vrtList,*state);
2417 }
2418 vrtList2.push_back(vID2);
2419 Trk::VertexID vID3;
2420 if(m_JXSubVtx) {
2421 // JX vertex
2422 if (m_constrJX && m_jxDaug_num>2) {
2423 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state,m_massJX);
2424 } else {
2425 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state);
2426 }
2427 vrtList2.push_back(vID3);
2428 // Mother vertex includes two subvertices: DisV and JX
2429 std::vector<const xAOD::TrackParticle*> tp;
2430 std::vector<double> tp_masses;
2431 if(m_constrMainV) {
2432 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2,*state,m_massMainV);
2433 } else {
2434 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2,*state);
2435 }
2436 }
2437 else { // m_JXSubVtx=false
2438 // Mother vertex includes just one subvertex (DisV) and JX tracks
2439 if(m_constrMainV) {
2440 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList2,*state,m_massMainV);
2441 } else {
2442 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList2,*state);
2443 }
2444 if (m_constrJX && m_jxDaug_num>2) {
2445 std::vector<Trk::VertexID> cnstV;
2446 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
2447 ATH_MSG_WARNING("addMassConstraint for JX failed");
2448 }
2449 }
2450 }
2451 if (m_constrJpsi) {
2452 std::vector<Trk::VertexID> cnstV;
2453 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
2454 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
2455 }
2456 }
2457 if (m_constrX && m_jxDaug_num==4) {
2458 std::vector<Trk::VertexID> cnstV;
2459 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksX,cnstV,*state,m_massX).isSuccess() ) {
2460 ATH_MSG_WARNING("addMassConstraint for X failed");
2461 }
2462 }
2463 // Do the work
2464 std::unique_ptr<Trk::VxCascadeInfo> fit_result = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state, pv_AOD.get(), pv_AOD.get() && m_firstDecayAtPV ? true : false) );
2465
2466 if (fit_result) {
2467 for(auto& v : fit_result->vertices()) {
2468 if(v->nTrackParticles()==0) {
2469 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
2470 v->setTrackParticleLinks(nullLinkVector);
2471 }
2472 }
2473 // reset links to original tracks
2474 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
2475
2476 // necessary to prevent memory leak
2477 fit_result->setSVOwnership(true);
2478
2479 // Chi2/DOF cut
2480 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
2481 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
2482
2483 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
2484 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
2485 size_t iMoth = cascadeVertices.size()-1;
2486 double lxy_SV1_sub = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
2487 double lxy_SV1 = m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[iMoth]);
2488
2489 if(chi2CutPassed && lxy_SV1>m_lxyDisV_cut && lxy_SV1_sub>m_lxyV0_cut) {
2490 chi2_V0_decor(*cascadeVertices[0]) = disVtx.V0vtx->chiSquared();
2491 ndof_V0_decor(*cascadeVertices[0]) = disVtx.V0vtx->numberDoF();
2492 if(disVtx.V0type==LAMBDA) {
2493 type_V0_decor(*cascadeVertices[0]) = "Lambda";
2494 }
2495 else if(disVtx.V0type==LAMBDABAR) {
2496 type_V0_decor(*cascadeVertices[0]) = "Lambdabar";
2497 }
2498 else if(disVtx.V0type==KS) {
2499 type_V0_decor(*cascadeVertices[0]) = "Ks";
2500 }
2501 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*disVtx.V0vtx) ? mAcc_gfit(*disVtx.V0vtx) : 0;
2502 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*disVtx.V0vtx) ? mAcc_gmass(*disVtx.V0vtx) : -1;
2503 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*disVtx.V0vtx) ? mAcc_gmasserr(*disVtx.V0vtx) : -1;
2504 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*disVtx.V0vtx) ? mAcc_gchisq(*disVtx.V0vtx) : 999999;
2505 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*disVtx.V0vtx) ? mAcc_gndof(*disVtx.V0vtx) : 0;
2506 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*disVtx.V0vtx) ? mAcc_gprob(*disVtx.V0vtx) : -1;
2507 trk_px.clear(); trk_py.clear(); trk_pz.clear();
2508 trk_px.push_back( disVtx.p4_V0track1.Px() ); trk_px.push_back( disVtx.p4_V0track2.Px() );
2509 trk_py.push_back( disVtx.p4_V0track1.Py() ); trk_py.push_back( disVtx.p4_V0track2.Py() );
2510 trk_pz.push_back( disVtx.p4_V0track1.Pz() ); trk_pz.push_back( disVtx.p4_V0track2.Pz() );
2511 trk_pxDeco(*cascadeVertices[0]) = trk_px;
2512 trk_pyDeco(*cascadeVertices[0]) = trk_py;
2513 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
2514 trk_px_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Px();
2515 trk_py_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Py();
2516 trk_pz_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Pz();
2517
2518 result.push_back( std::make_pair(fit_result.release(),nullptr) );
2519 }
2520 }
2521 } // m_extraTrk1MassHypo<=0
2522 else { // m_extraTrk1MassHypo>0
2523 std::vector<double> massesExtra{m_extraTrk1MassHypo};
2524 std::vector<double> massesJXExtra = massesJX; massesJXExtra.push_back(m_extraTrk1MassHypo);
2525
2526 for(const xAOD::TrackParticle* tpExtra : *trackContainer) {
2527 if ( tpExtra->pt()<m_extraTrk1MinPt ) continue;
2528 if ( !m_trkSelector->decision(*tpExtra, nullptr) ) continue;
2529 // Check identical tracks in input
2530 if(std::find(tracksJX.cbegin(),tracksJX.cend(),tpExtra) != tracksJX.cend()) continue;
2531 if(std::find(tracksV0.cbegin(),tracksV0.cend(),tpExtra) != tracksV0.cend()) continue;
2532 if(tpExtra == disVtx.track) continue;
2533
2534 TLorentzVector tmp;
2535 tmp.SetPtEtaPhiM(tpExtra->pt(),tpExtra->eta(),tpExtra->phi(),m_extraTrk1MassHypo);
2536 double main_mass = (p4_moth+tmp).M();
2537 if(m_useImprovedMass) {
2538 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_disV).M() + m_massJpsi;
2539 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_disV).M() + m_massJX;
2540 if(m_massDisV>0) main_mass += - p4_disV.M() + m_massDisV;
2541 }
2542 if (main_mass < m_MassLower || main_mass > m_MassUpper) continue;
2543
2544 std::vector<const xAOD::TrackParticle*> tracksExtra{tpExtra};
2545 std::vector<const xAOD::TrackParticle*> tracksJXExtra = tracksJX; tracksJXExtra.push_back(tpExtra);
2546
2547 // Apply the user's settings to the fitter
2548 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
2549 // Robustness: http://cdsweb.cern.ch/record/685551
2550 int robustness = 0;
2551 m_iVertexFitter->setRobustness(robustness, *state);
2552 // Build up the topology
2553 // Vertex list
2554 std::vector<Trk::VertexID> vrtList;
2555 std::vector<Trk::VertexID> vrtList2;
2556 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
2557 // V0 vertex
2558 Trk::VertexID vID1;
2559 if (m_constrV0) {
2560 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,disVtx.V0type==KS?m_massKs:m_massLd);
2561 } else {
2562 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
2563 }
2564 vrtList.push_back(vID1);
2565 // Displaced vertex
2566 Trk::VertexID vID2;
2567 if (m_constrDisV) {
2568 vID2 = m_iVertexFitter->nextVertex(tracks3,massesDis3,vrtList,*state,m_massDisV);
2569 } else {
2570 vID2 = m_iVertexFitter->nextVertex(tracks3,massesDis3,vrtList,*state);
2571 }
2572 vrtList2.push_back(vID2);
2573 Trk::VertexID vID3;
2574 if(m_JXSubVtx) {
2575 // JXExtra vertex
2576 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state);
2577 vrtList2.push_back(vID3);
2578 // Mother vertex
2579 if(m_constrMainV) {
2580 m_iVertexFitter->nextVertex(tracksExtra,massesExtra,vrtList2,*state,m_massMainV);
2581 } else {
2582 m_iVertexFitter->nextVertex(tracksExtra,massesExtra,vrtList2,*state);
2583 }
2584 }
2585 else { // m_JXSubVtx=false
2586 // Mother vertex includes just one subvertex (DisV) and JX tracks + extra track
2587 if(m_constrMainV) {
2588 vID3 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList2,*state,m_massMainV);
2589 } else {
2590 vID3 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList2,*state);
2591 }
2592 }
2593 if (m_constrJX && m_jxDaug_num>2) {
2594 std::vector<Trk::VertexID> cnstV;
2595 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
2596 ATH_MSG_WARNING("addMassConstraint for JX failed");
2597 }
2598 }
2599 if (m_constrJpsi) {
2600 std::vector<Trk::VertexID> cnstV;
2601 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
2602 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
2603 }
2604 }
2605 if (m_constrX && m_jxDaug_num==4) {
2606 std::vector<Trk::VertexID> cnstV;
2607 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksX,cnstV,*state,m_massX).isSuccess() ) {
2608 ATH_MSG_WARNING("addMassConstraint for X failed");
2609 }
2610 }
2611 // Do the work
2612 std::unique_ptr<Trk::VxCascadeInfo> fit_result = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state, pv_AOD.get(), pv_AOD.get() && m_firstDecayAtPV ? true : false) );
2613
2614 if (fit_result) {
2615 for(auto& v : fit_result->vertices()) {
2616 if(v->nTrackParticles()==0) {
2617 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
2618 v->setTrackParticleLinks(nullLinkVector);
2619 }
2620 }
2621 // reset links to original tracks
2622 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
2623
2624 // necessary to prevent memory leak
2625 fit_result->setSVOwnership(true);
2626
2627 // Chi2/DOF cut
2628 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
2629 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
2630
2631 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
2632 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
2633 size_t iMoth = cascadeVertices.size()-1;
2634 double lxy_SV1_sub = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
2635 double lxy_SV1 = m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[iMoth]);
2636
2637 if(chi2CutPassed && lxy_SV1>m_lxyDisV_cut && lxy_SV1_sub>m_lxyV0_cut) {
2638 chi2_V0_decor(*cascadeVertices[0]) = disVtx.V0vtx->chiSquared();
2639 ndof_V0_decor(*cascadeVertices[0]) = disVtx.V0vtx->numberDoF();
2640 if(disVtx.V0type==LAMBDA) {
2641 type_V0_decor(*cascadeVertices[0]) = "Lambda";
2642 }
2643 else if(disVtx.V0type==LAMBDABAR) {
2644 type_V0_decor(*cascadeVertices[0]) = "Lambdabar";
2645 }
2646 else if(disVtx.V0type==KS) {
2647 type_V0_decor(*cascadeVertices[0]) = "Ks";
2648 }
2649 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*disVtx.V0vtx) ? mAcc_gfit(*disVtx.V0vtx) : 0;
2650 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*disVtx.V0vtx) ? mAcc_gmass(*disVtx.V0vtx) : -1;
2651 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*disVtx.V0vtx) ? mAcc_gmasserr(*disVtx.V0vtx) : -1;
2652 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*disVtx.V0vtx) ? mAcc_gchisq(*disVtx.V0vtx) : 999999;
2653 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*disVtx.V0vtx) ? mAcc_gndof(*disVtx.V0vtx) : 0;
2654 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*disVtx.V0vtx) ? mAcc_gprob(*disVtx.V0vtx) : -1;
2655 trk_px.clear(); trk_py.clear(); trk_pz.clear();
2656 trk_px.push_back( disVtx.p4_V0track1.Px() ); trk_px.push_back( disVtx.p4_V0track2.Px() );
2657 trk_py.push_back( disVtx.p4_V0track1.Py() ); trk_py.push_back( disVtx.p4_V0track2.Py() );
2658 trk_pz.push_back( disVtx.p4_V0track1.Pz() ); trk_pz.push_back( disVtx.p4_V0track2.Pz() );
2659 trk_pxDeco(*cascadeVertices[0]) = trk_px;
2660 trk_pyDeco(*cascadeVertices[0]) = trk_py;
2661 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
2662 trk_px_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Px();
2663 trk_py_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Py();
2664 trk_pz_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Pz();
2665
2666 result.push_back( std::make_pair(fit_result.release(),nullptr) );
2667 }
2668 }
2669 } // loop over trackContainer
2670 } // m_extraTrk1MassHypo>0
2671
2672 if(pv_xAOD) {
2673 for(auto cascade_info_pair : result) {
2674 if(cascade_info_pair.first && cascade_info_pair.first->getParticleMoms().size()>0) {
2675 size_t index = cascade_info_pair.first->getParticleMoms().size() - 1;
2676 const std::vector<TLorentzVector> &mom = cascade_info_pair.first->getParticleMoms()[index];
2677 const Amg::MatrixX &cov = cascade_info_pair.first->getCovariance()[index];
2678 const xAOD::Vertex* mainVertex = cascade_info_pair.first->vertices()[index];
2679 xAOD::BPhysHypoHelper vtx(m_hypoName, mainVertex);
2680 bool isInDefaultPVCont = false;
2681 for(const xAOD::Vertex* pvVtx : *defaultPVContainer) {
2682 if(pv_xAOD == pvVtx) { isInDefaultPVCont = true; break; }
2683 }
2684 if(isInDefaultPVCont) vtx.setPv( pv_xAOD, defaultPVContainer, pvtype );
2685 else vtx.setPv( pv_xAOD, pvContainer, pvtype );
2686 if(origPv_xAOD) vtx.setOrigPv( origPv_xAOD, defaultPVContainer, pvtype );
2687 vtx.setLxy ( m_CascadeTools->lxy (mom, vtx.vtx(), pv_xAOD), pvtype );
2688 vtx.setLxyErr ( m_CascadeTools->lxyError (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2689 vtx.setA0 ( m_CascadeTools->a0 (mom, vtx.vtx(), pv_xAOD), pvtype );
2690 vtx.setA0Err ( m_CascadeTools->a0Error (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2691 vtx.setA0xy ( m_CascadeTools->a0xy (mom, vtx.vtx(), pv_xAOD), pvtype );
2692 vtx.setA0xyErr( m_CascadeTools->a0xyError(mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2693 vtx.setZ0 ( m_CascadeTools->a0z (mom, vtx.vtx(), pv_xAOD), pvtype );
2694 vtx.setZ0Err ( m_CascadeTools->a0zError (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2695 vtx.setRefitPVStatus( 0, pvtype );
2696 // Proper decay times
2697 vtx.setTau( m_CascadeTools->tau(mom, vtx.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
2698 vtx.setTauErr( m_CascadeTools->tauError(mom, cov, vtx.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
2699 vtx.setTau( m_CascadeTools->tau(mom, vtx.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
2700 vtx.setTauErr( m_CascadeTools->tauError(mom, cov, vtx.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
2701 }
2702 }
2703 }
2704
2705 return result;
2706 }
2707
2708 void JpsiXPlusDisplaced::fitV0Container(const EventContext& ctx, xAOD::VertexContainer* V0ContainerNew, const std::vector<const xAOD::TrackParticle*>& selectedTracks, const std::vector<const xAOD::TrackParticleContainer*>& trackCols) const {
2709
2710 SG::AuxElement::Decorator<std::string> mDec_type("Type_V0Vtx");
2711 SG::AuxElement::Decorator<int> mDec_gfit("gamma_fit");
2712 SG::AuxElement::Decorator<float> mDec_gmass("gamma_mass");
2713 SG::AuxElement::Decorator<float> mDec_gmasserr("gamma_massError");
2714 SG::AuxElement::Decorator<float> mDec_gchisq("gamma_chisq");
2715 SG::AuxElement::Decorator<int> mDec_gndof("gamma_ndof");
2716 SG::AuxElement::Decorator<float> mDec_gprob("gamma_probability");
2717
2718 std::vector<const xAOD::TrackParticle*> posTracks;
2719 std::vector<const xAOD::TrackParticle*> negTracks;
2720 for(const xAOD::TrackParticle* TP : selectedTracks) {
2721 if(TP->charge()>0) posTracks.push_back(TP);
2722 else negTracks.push_back(TP);
2723 }
2724
2725 for(const xAOD::TrackParticle* TP1 : posTracks) {
2726 const Trk::Perigee& aPerigee1 = TP1->perigeeParameters();
2727 for(const xAOD::TrackParticle* TP2 : negTracks) {
2728 const Trk::Perigee& aPerigee2 = TP2->perigeeParameters();
2729 int sflag(0), errorcode(0);
2730 Amg::Vector3D startingPoint = m_vertexEstimator->getCirclesIntersectionPoint(&aPerigee1,&aPerigee2,sflag,errorcode);
2731 if (errorcode != 0) {startingPoint(0) = 0.0; startingPoint(1) = 0.0; startingPoint(2) = 0.0;}
2732
2733 if (errorcode == 0 || errorcode == 5 || errorcode == 6 || errorcode == 8) {
2734 Trk::PerigeeSurface perigeeSurface(startingPoint);
2735 const Trk::TrackParameters* extrapolatedPerigee1 = m_extrapolator->extrapolate(ctx,TP1->perigeeParameters(), perigeeSurface).release();
2736 const Trk::TrackParameters* extrapolatedPerigee2 = m_extrapolator->extrapolate(ctx,TP2->perigeeParameters(), perigeeSurface).release();
2737 std::vector<std::unique_ptr<const Trk::TrackParameters> > cleanup;
2738 if(!extrapolatedPerigee1) extrapolatedPerigee1 = &TP1->perigeeParameters();
2739 else cleanup.push_back(std::unique_ptr<const Trk::TrackParameters>(extrapolatedPerigee1));
2740 if(!extrapolatedPerigee2) extrapolatedPerigee2 = &TP2->perigeeParameters();
2741 else cleanup.push_back(std::unique_ptr<const Trk::TrackParameters>(extrapolatedPerigee2));
2742 if(extrapolatedPerigee1 && extrapolatedPerigee2) {
2743 bool pass = false;
2744 TLorentzVector v1; TLorentzVector v2;
2745 if(!pass) {
2746 v1.SetXYZM(extrapolatedPerigee1->momentum().x(),extrapolatedPerigee1->momentum().y(),extrapolatedPerigee1->momentum().z(),m_mass_proton);
2747 v2.SetXYZM(extrapolatedPerigee2->momentum().x(),extrapolatedPerigee2->momentum().y(),extrapolatedPerigee2->momentum().z(),m_mass_pion);
2748 if((v1+v2).M()>1030.0 && (v1+v2).M()<1200.0) pass = true;
2749 }
2750 if(!pass) {
2751 v1.SetXYZM(extrapolatedPerigee1->momentum().x(),extrapolatedPerigee1->momentum().y(),extrapolatedPerigee1->momentum().z(),m_mass_pion);
2752 v2.SetXYZM(extrapolatedPerigee2->momentum().x(),extrapolatedPerigee2->momentum().y(),extrapolatedPerigee2->momentum().z(),m_mass_proton);
2753 if((v1+v2).M()>1030.0 && (v1+v2).M()<1200.0) pass = true;
2754 }
2755 if(!pass) {
2756 v1.SetXYZM(extrapolatedPerigee1->momentum().x(),extrapolatedPerigee1->momentum().y(),extrapolatedPerigee1->momentum().z(),m_mass_pion);
2757 v2.SetXYZM(extrapolatedPerigee2->momentum().x(),extrapolatedPerigee2->momentum().y(),extrapolatedPerigee2->momentum().z(),m_mass_pion);
2758 if((v1+v2).M()>430.0 && (v1+v2).M()<565.0) pass = true;
2759 }
2760 if(pass) {
2761 std::vector<const xAOD::TrackParticle*> tracksV0;
2762 tracksV0.push_back(TP1); tracksV0.push_back(TP2);
2763 std::unique_ptr<xAOD::Vertex> V0vtx = m_iV0Fitter->fit(ctx, tracksV0, startingPoint);
2764 if(V0vtx && V0vtx->chiSquared()>=0) {
2765 double chi2DOF = V0vtx->chiSquared()/V0vtx->numberDoF();
2766 if(chi2DOF>m_chi2cut_V0) continue;
2767
2768 double massSig_V0_Lambda1 = std::abs(m_V0Tools->invariantMass(V0vtx.get(), m_massesV0_ppi)-m_mass_Lambda)/m_V0Tools->invariantMassError(V0vtx.get(), m_massesV0_ppi);
2769 double massSig_V0_Lambda2 = std::abs(m_V0Tools->invariantMass(V0vtx.get(), m_massesV0_pip)-m_mass_Lambda)/m_V0Tools->invariantMassError(V0vtx.get(), m_massesV0_pip);
2770 double massSig_V0_Ks = std::abs(m_V0Tools->invariantMass(V0vtx.get(), m_massesV0_pipi)-m_mass_Ks)/m_V0Tools->invariantMassError(V0vtx.get(), m_massesV0_pipi);
2771 if(massSig_V0_Lambda1<=massSig_V0_Lambda2 && massSig_V0_Lambda1<=massSig_V0_Ks) {
2772 mDec_type(*V0vtx.get()) = "Lambda";
2773 }
2774 else if(massSig_V0_Lambda2<=massSig_V0_Lambda1 && massSig_V0_Lambda2<=massSig_V0_Ks) {
2775 mDec_type(*V0vtx.get()) = "Lambdabar";
2776 }
2777 else if(massSig_V0_Ks<=massSig_V0_Lambda1 && massSig_V0_Ks<=massSig_V0_Lambda2) {
2778 mDec_type(*V0vtx.get()) = "Ks";
2779 }
2780
2781 int gamma_fit = 0; int gamma_ndof = 0; double gamma_chisq = 999999.;
2782 double gamma_prob = -1., gamma_mass = -1., gamma_massErr = -1.;
2783 std::unique_ptr<xAOD::Vertex> gammaVtx = m_iGammaFitter->fit(ctx, tracksV0, m_V0Tools->vtx(V0vtx.get()));
2784 if (gammaVtx) {
2785 gamma_fit = 1;
2786 gamma_mass = m_V0Tools->invariantMass(gammaVtx.get(),m_mass_e,m_mass_e);
2787 gamma_massErr = m_V0Tools->invariantMassError(gammaVtx.get(),m_mass_e,m_mass_e);
2788 gamma_chisq = m_V0Tools->chisq(gammaVtx.get());
2789 gamma_ndof = m_V0Tools->ndof(gammaVtx.get());
2790 gamma_prob = m_V0Tools->vertexProbability(gammaVtx.get());
2791 }
2792 mDec_gfit(*V0vtx.get()) = gamma_fit;
2793 mDec_gmass(*V0vtx.get()) = gamma_mass;
2794 mDec_gmasserr(*V0vtx.get()) = gamma_massErr;
2795 mDec_gchisq(*V0vtx.get()) = gamma_chisq;
2796 mDec_gndof(*V0vtx.get()) = gamma_ndof;
2797 mDec_gprob(*V0vtx.get()) = gamma_prob;
2798
2799 xAOD::BPhysHelper V0_helper(V0vtx.get());
2800 V0_helper.setRefTrks(); // AOD only method
2801
2802 if(not trackCols.empty()){
2803 try {
2804 JpsiUpsilonCommon::RelinkVertexTracks(trackCols, V0vtx.get());
2805 } catch (std::runtime_error const& e) {
2806 ATH_MSG_ERROR(e.what());
2807 return;
2808 }
2809 }
2810
2811 V0ContainerNew->push_back(std::move(V0vtx));
2812 }
2813 }
2814 }
2815 }
2816 }
2817 }
2818 }
2819
2820 template<size_t NTracks>
2821 const xAOD::Vertex* JpsiXPlusDisplaced::FindVertex(const xAOD::VertexContainer* cont, const xAOD::Vertex* v) const {
2822 for (const xAOD::Vertex* v1 : *cont) {
2823 assert(v1->nTrackParticles() == NTracks);
2824 std::array<const xAOD::TrackParticle*, NTracks> a1;
2825 std::array<const xAOD::TrackParticle*, NTracks> a2;
2826 for(size_t i=0; i<NTracks; i++){
2827 a1[i] = v1->trackParticle(i);
2828 a2[i] = v->trackParticle(i);
2829 }
2830 std::sort(a1.begin(), a1.end());
2831 std::sort(a2.begin(), a2.end());
2832 if(a1 == a2) return v1;
2833 }
2834 return nullptr;
2835 }
2836}
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_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.
passed
bool passed(DecisionID id, const DecisionIDContainer &)
checks if required decision ID is in the set of IDs in the container
Definition TrigCompositeUtilsRoot.cxx:118
a
static Double_t a
Definition LArPhysWaveHECTool.cxx:38
size
size_t size() const
Number of registered mappings.
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
DerivationFramework::BPhysPVCascadeTools::PrepareVertexLinks
static void PrepareVertexLinks(Trk::VxCascadeInfo *result, const xAOD::TrackParticleContainer *importedTrackCollection)
Definition BPhysPVCascadeTools.cxx:204
DerivationFramework::BPhysPVCascadeTools::getParticleMass
static double getParticleMass(const HepPDT::ParticleDataTable *pdt, int pdg)
Definition BPhysPVCascadeTools.cxx:491
DerivationFramework::BPhysPVCascadeTools::SetVectorInfo
static void SetVectorInfo(xAOD::BPhysHelper &, const Trk::VxCascadeInfo *)
Definition BPhysPVCascadeTools.cxx:424
DerivationFramework::JpsiUpsilonCommon::RelinkVertexTracks
static void RelinkVertexTracks(std::span< const xAOD::TrackParticleContainer *const > trkcols, xAOD::Vertex *vtx)
Definition JpsiUpsilonCommon.cxx:129
DerivationFramework::JpsiXPlusDisplaced::MesonCandidateVector::vector
const std::vector< MesonCandidate > & vector() const
Definition JpsiXPlusDisplaced.cxx:52
DerivationFramework::JpsiXPlusDisplaced::m_trkSelector
ToolHandle< Trk::ITrackSelectorTool > m_trkSelector
Definition JpsiXPlusDisplaced.h:190
DerivationFramework::JpsiXPlusDisplaced::m_eventInfo_key
SG::ReadHandleKey< xAOD::EventInfo > m_eventInfo_key
Definition JpsiXPlusDisplaced.h:95
DerivationFramework::JpsiXPlusDisplaced::m_iV0Fitter
ToolHandle< Trk::TrkV0VertexFitter > m_iV0Fitter
Definition JpsiXPlusDisplaced.h:185
DerivationFramework::JpsiXPlusDisplaced::m_cascadeOutputKeys
SG::WriteHandleKeyArray< xAOD::VertexContainer > m_cascadeOutputKeys
Definition JpsiXPlusDisplaced.h:88
DerivationFramework::JpsiXPlusDisplaced::getD0Candidate
MesonCandidate getD0Candidate(const EventContext &ctx, const xAOD::Vertex *JXvtx, const xAOD::TrackParticle *extraTrk1, const xAOD::TrackParticle *extraTrk2) const
Definition JpsiXPlusDisplaced.cxx:1281
DerivationFramework::JpsiXPlusDisplaced::m_TrkParticleCollection
SG::ReadHandleKey< xAOD::TrackParticleContainer > m_TrkParticleCollection
Definition JpsiXPlusDisplaced.h:91
DerivationFramework::JpsiXPlusDisplaced::m_iVertexFitter
ToolHandle< Trk::TrkVKalVrtFitter > m_iVertexFitter
Definition JpsiXPlusDisplaced.h:184
DerivationFramework::JpsiXPlusDisplaced::m_iGammaFitter
ToolHandle< Trk::IVertexFitter > m_iGammaFitter
Definition JpsiXPlusDisplaced.h:186
DerivationFramework::JpsiXPlusDisplaced::m_CascadeTools
PublicToolHandle< DerivationFramework::CascadeTools > m_CascadeTools
Definition JpsiXPlusDisplaced.h:192
DerivationFramework::JpsiXPlusDisplaced::JpsiXPlusDisplaced
JpsiXPlusDisplaced(const std::string &type, const std::string &name, const IInterface *parent)
Definition JpsiXPlusDisplaced.cxx:56
DerivationFramework::JpsiXPlusDisplaced::m_vertexJXHypoNames
std::vector< std::string > m_vertexJXHypoNames
Definition JpsiXPlusDisplaced.h:87
DerivationFramework::JpsiXPlusDisplaced::fitMainVtx
std::vector< std::pair< Trk::VxCascadeInfo *, Trk::VxCascadeInfo * > > fitMainVtx(const EventContext &ctx, const xAOD::Vertex *JXvtx, const std::vector< double > &massesJX, const xAOD::Vertex *V0vtx, const V0Enum V0, const xAOD::TrackParticleContainer *trackContainer, const std::vector< const xAOD::TrackParticleContainer * > &trackCols, const xAOD::VertexContainer *defaultPVContainer, const xAOD::VertexContainer *pvContainer) const
Definition JpsiXPlusDisplaced.cxx:1312
DerivationFramework::JpsiXPlusDisplaced::m_cascadeOutputKeys_mvc
SG::WriteHandleKeyArray< xAOD::VertexContainer > m_cascadeOutputKeys_mvc
Definition JpsiXPlusDisplaced.h:89
DerivationFramework::JpsiXPlusDisplaced::m_RelinkContainers
SG::ReadHandleKeyArray< xAOD::TrackParticleContainer > m_RelinkContainers
Definition JpsiXPlusDisplaced.h:96
DerivationFramework::JpsiXPlusDisplaced::m_refPVContainerName
SG::WriteHandleKey< xAOD::VertexContainer > m_refPVContainerName
Definition JpsiXPlusDisplaced.h:94
DerivationFramework::JpsiXPlusDisplaced::m_vertexV0ContainerKey
SG::ReadHandleKey< xAOD::VertexContainer > m_vertexV0ContainerKey
Definition JpsiXPlusDisplaced.h:86
DerivationFramework::JpsiXPlusDisplaced::getDpmCandidate
MesonCandidate getDpmCandidate(const EventContext &ctx, const xAOD::Vertex *JXvtx, const xAOD::TrackParticle *extraTrk1, const xAOD::TrackParticle *extraTrk2, const xAOD::TrackParticle *extraTrk3) const
Definition JpsiXPlusDisplaced.cxx:1245
DerivationFramework::JpsiXPlusDisplaced::FindVertex
const xAOD::Vertex * FindVertex(const xAOD::VertexContainer *cont, const xAOD::Vertex *v) const
Definition JpsiXPlusDisplaced.cxx:2821
DerivationFramework::JpsiXPlusDisplaced::getXiCandidate
XiCandidate getXiCandidate(const EventContext &ctx, const xAOD::Vertex *V0vtx, const V0Enum V0, const xAOD::TrackParticle *track3) const
Definition JpsiXPlusDisplaced.cxx:1178
DerivationFramework::JpsiXPlusDisplaced::m_v0VtxOutputKey
SG::WriteHandleKey< xAOD::VertexContainer > m_v0VtxOutputKey
Definition JpsiXPlusDisplaced.h:90
DerivationFramework::JpsiXPlusDisplaced::fitV0Container
void fitV0Container(const EventContext &ctx, xAOD::VertexContainer *V0ContainerNew, const std::vector< const xAOD::TrackParticle * > &selectedTracks, const std::vector< const xAOD::TrackParticleContainer * > &trackCols) const
Definition JpsiXPlusDisplaced.cxx:2708
DerivationFramework::JpsiXPlusDisplaced::m_vertexJXContainerKey
SG::ReadHandleKey< xAOD::VertexContainer > m_vertexJXContainerKey
Definition JpsiXPlusDisplaced.h:85
DerivationFramework::JpsiXPlusDisplaced::addBranches
virtual StatusCode addBranches(const EventContext &ctx) const override
Definition JpsiXPlusDisplaced.cxx:549
DerivationFramework::JpsiXPlusDisplaced::m_VxPrimaryCandidateName
SG::ReadHandleKey< xAOD::VertexContainer > m_VxPrimaryCandidateName
Definition JpsiXPlusDisplaced.h:92
DerivationFramework::JpsiXPlusDisplaced::m_v0TrkSelector
ToolHandle< Trk::ITrackSelectorTool > m_v0TrkSelector
Definition JpsiXPlusDisplaced.h:191
DerivationFramework::JpsiXPlusDisplaced::m_pvRefitter
PublicToolHandle< Analysis::PrimaryVertexRefitter > m_pvRefitter
Definition JpsiXPlusDisplaced.h:187
DerivationFramework::JpsiXPlusDisplaced::d0Pass
bool d0Pass(const EventContext &ctx, const xAOD::TrackParticle *track, const xAOD::Vertex *PV) const
Definition JpsiXPlusDisplaced.cxx:1168
DerivationFramework::JpsiXPlusDisplaced::m_pvContainerName
SG::ReadHandleKey< xAOD::VertexContainer > m_pvContainerName
Definition JpsiXPlusDisplaced.h:93
DerivationFramework::JpsiXPlusDisplaced::performSearch
StatusCode performSearch(std::vector< std::pair< Trk::VxCascadeInfo *, Trk::VxCascadeInfo * > > &cascadeinfoContainer, const std::vector< std::pair< const xAOD::Vertex *, V0Enum > > &selectedV0Candidates, const std::vector< const xAOD::TrackParticle * > &tracksDisplaced, const EventContext &ctx) const
Definition JpsiXPlusDisplaced.cxx:384
DerivationFramework::JpsiXPlusDisplaced::m_trackToVertexTool
ToolHandle< Reco::ITrackToVertex > m_trackToVertexTool
Definition JpsiXPlusDisplaced.h:189
DerivationFramework::JpsiXPlusDisplaced::m_extrapolator
ToolHandle< Trk::IExtrapolator > m_extrapolator
Definition JpsiXPlusDisplaced.h:194
DerivationFramework::JpsiXPlusDisplaced::m_partPropSvc
ServiceHandle< IPartPropSvc > m_partPropSvc
Definition JpsiXPlusDisplaced.h:182
DerivationFramework::JpsiXPlusDisplaced::m_vertexEstimator
ToolHandle< InDet::VertexPointEstimator > m_vertexEstimator
Definition JpsiXPlusDisplaced.h:193
DerivationFramework::JpsiXPlusDisplaced::fitTracks
std::unique_ptr< xAOD::Vertex > fitTracks(const EventContext &ctx, const xAOD::TrackParticle *track1, const xAOD::TrackParticle *track2, const xAOD::TrackParticle *track3=nullptr) const
Definition JpsiXPlusDisplaced.cxx:1228
DerivationFramework::JpsiXPlusDisplaced::m_V0Tools
PublicToolHandle< Trk::V0Tools > m_V0Tools
Definition JpsiXPlusDisplaced.h:188
SG::ReadHandleKey
Property holding a SG store/key/clid from which a ReadHandle is made.
Definition StoreGate/StoreGate/ReadHandleKey.h:40
SG::ReadHandle::ptr
const_pointer_type ptr()
Dereference the pointer.
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::WriteHandleKey
Property holding a SG store/key/clid from which a WriteHandle is made.
Definition StoreGate/StoreGate/WriteHandleKey.h:40
SG::WriteHandle::cptr
const_pointer_type cptr() const
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.
Trk::ParametersBase::momentum
const Amg::Vector3D & momentum() const
Access method for the momentum.
Trk::PerigeeSurface
Class describing the Line to which the Perigee refers to.
Definition PerigeeSurface.h:43
xAOD::BPhysHelper::setZ0
float setZ0(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
longitudinal impact parameter
Definition BPhysHelper.cxx:1067
xAOD::BPhysHelper::origPv
const xAOD::Vertex * origPv(const pv_type vertexType=BPhysHelper::PV_MIN_A0)
original PV
Definition BPhysHelper.cxx:807
xAOD::BPhysHelper::setLxyErr
bool setLxyErr(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
its error
Definition BPhysHelper.cxx:899
xAOD::BPhysHelper::setLxy
bool setLxy(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
Set the transverse decay distance and its error measured between the refitted primary vertex of type ...
Definition BPhysHelper.cxx:888
xAOD::BPhysHelper::setZ0Err
float setZ0Err(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
longitudinal impact parameter error
Definition BPhysHelper.cxx:1078
xAOD::BPhysHelper::vtx
const xAOD::Vertex * vtx() const
Getter method for the cached xAOD::Vertex.
Definition BPhysHelper.h:108
xAOD::BPhysHelper::setA0
float setA0(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
Set the 3D and transverse impact parameters and their error.
Definition BPhysHelper.cxx:1023
xAOD::BPhysHelper::refTrks
const std::vector< TVector3 > & refTrks()
Returns refitted track momenta.
Definition BPhysHelper.cxx:142
xAOD::BPhysHelper::setPv
bool setPv(const xAOD::Vertex *pv, const xAOD::VertexContainer *vertexContainer, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
Set the refitted collision vertex of type pv_type.
Definition BPhysHelper.cxx:831
xAOD::BPhysHelper::pv
const xAOD::Vertex * pv(const pv_type vertexType=BPhysHelper::PV_MIN_A0)
Get the refitted collision vertex of type pv_type.
Definition BPhysHelper.cxx:796
xAOD::BPhysHelper::refTrk
TVector3 refTrk(const size_t index)
Returns i-th refitted track 3-momentum.
Definition BPhysHelper.cxx:126
xAOD::BPhysHelper::pv_type
pv_type
: Enum type of the PV
Definition BPhysHelper.h:475
xAOD::BPhysHelper::setRefTrks
bool setRefTrks(std::vector< float > px, std::vector< float > py, std::vector< float > pz)
Sets refitted track momenta.
Definition BPhysHelper.cxx:286
xAOD::BPhysHelper::setA0Err
float setA0Err(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
3D impact parameter error
Definition BPhysHelper.cxx:1034
xAOD::BPhysHelper::setRefitPVStatus
bool setRefitPVStatus(int code, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
Set the exitCode of the refitter for vertex of type pv_type.
Definition BPhysHelper.cxx:844
xAOD::BPhysHelper::setA0xyErr
float setA0xyErr(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
transverse impact parameter error
Definition BPhysHelper.cxx:1056
xAOD::BPhysHelper::nRefTrks
int nRefTrks()
Returns number of stored refitted track momenta.
Definition BPhysHelper.cxx:115
xAOD::BPhysHelper::setA0xy
float setA0xy(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
transverse impact parameter
Definition BPhysHelper.cxx:1045
xAOD::BPhysHelper::setOrigPv
bool setOrigPv(const xAOD::Vertex *pv, const xAOD::VertexContainer *vertexContainer, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
Set the original collision vertex of type pv_type.
Definition BPhysHelper.cxx:818
xAOD::BPhysHypoHelper::setMass
bool setMass(const float val)
Set given invariant mass and its error.
Definition BPhysHypoHelper.cxx:52
xAOD::BPhysHypoHelper::setPass
bool setPass(bool passVal)
get the pass flag for this hypothesis
Definition BPhysHypoHelper.cxx:367
xAOD::BPhysHypoHelper::setTau
bool setTau(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0, const tau_type tauType=BPhysHypoHelper::TAU_CONST_MASS)
: Set the proper decay time and error.
Definition BPhysHypoHelper.cxx:143
xAOD::BPhysHypoHelper::setMassErr
bool setMassErr(const float val)
invariant mass error
Definition BPhysHypoHelper.cxx:57
xAOD::BPhysHypoHelper::setTauErr
bool setTauErr(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0, const tau_type tauType=BPhysHypoHelper::TAU_CONST_MASS)
proper decay time error
Definition BPhysHypoHelper.cxx:179
xAOD::TrackParticle_v1::perigeeParameters
const Trk::Perigee & perigeeParameters() const
Returns the Trk::MeasuredPerigee track parameters.
Definition TrackParticle_v1.cxx:437
xAOD::TrackParticle_v1::phi
virtual double phi() const override final
The azimuthal angle ( ) of the particle (has range to .).
xAOD::TrackParticle_v1::pt
virtual double pt() const override final
The transverse momentum ( ) of the particle.
Definition TrackParticle_v1.cxx:75
xAOD::TrackParticle_v1::eta
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
Definition TrackParticle_v1.cxx:79
xAOD::Vertex_v1::nTrackParticles
size_t nTrackParticles() const
Get the number of tracks associated with this vertex.
Definition Vertex_v1.cxx:292
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::numberDoF
float numberDoF() const
Returns the number of degrees of freedom of the vertex fit as float.
xAOD::Vertex_v1::chiSquared
float chiSquared() const
Returns the of the vertex fit as float.
xAOD::Vertex_v1::position
const Amg::Vector3D & position() const
Returns the 3-pos.
Definition Vertex_v1.cxx:106
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition EventPrimitives.h:27
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47
DerivationFramework
THE reconstruction tool.
Definition ParticleSortingAlg.h:24
DerivationFramework::VertexLink
ElementLink< xAOD::VertexContainer > VertexLink
Definition JpsiPlusDpstCascade.cxx:23
DerivationFramework::VertexLinkVector
std::vector< VertexLink > VertexLinkVector
Definition JpsiPlusDpstCascade.cxx:24
MC::PSI2S
static const int PSI2S
Definition HepMCHelpers.h:130
MC::MUON
static const int MUON
Definition HepMCHelpers.h:78
MC::DPLUS
static const int DPLUS
Definition HepMCHelpers.h:118
MC::BCPLUS
static const int BCPLUS
Definition HepMCHelpers.h:124
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::D0
static const int D0
Definition HepMCHelpers.h:120
MC::PROTON
static const int PROTON
Definition HepMCHelpers.h:125
Trk::Perigee
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
Definition Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:33
Trk::pz
@ pz
global momentum (cartesian)
Definition ParamDefs.h:61
Trk::d0
@ d0
Definition ParamDefs.h:63
Trk::px
@ px
Definition ParamDefs.h:59
Trk::py
@ py
Definition ParamDefs.h:60
Trk::TrackParameters
ParametersBase< TrackParametersDim, Charged > TrackParameters
Definition Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:27
get_generator_info.result
result
Definition get_generator_info.py:21
index
Definition index.py:1
std::sort
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
Definition DVL_algorithms.h:554
xAOD::TrackParticle
TrackParticle_v1 TrackParticle
Reference the current persistent version:
Definition Event/xAOD/xAODTracking/xAODTracking/TrackParticle.h:13
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
xAOD::TrackParticleContainer
TrackParticleContainer_v1 TrackParticleContainer
Definition of the current "TrackParticle container version".
Definition Event/xAOD/xAODTracking/xAODTracking/TrackParticleContainer.h:14
xAOD::numberOfSCTHits
@ numberOfSCTHits
number of hits in SCT [unit8_t].
Definition TrackingPrimitives.h:269
xAOD::numberOfPixelHits
@ numberOfPixelHits
these are the pixel hits, including the b-layer [unit8_t].
Definition TrackingPrimitives.h:260
type
Generated on for ATLAS Offline Software by doxygen 1.16.1