ATLAS Offline Software
Loading...
Searching...
No Matches
JpsiXPlusDisplaced.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2025 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", this),
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(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(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(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(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(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 xAOD::TrackParticle* track, const xAOD::Vertex* PV) const {
1169 bool pass = false;
1170 const EventContext& ctx = Gaudi::Hive::currentContext();
1171 std::unique_ptr<Trk::Perigee> per = m_trackToVertexTool->perigeeAtVertex(ctx, *track, PV->position());
1172 if(!per) return pass;
1173 double d0 = per->parameters()[Trk::d0];
1174 double sig_d0 = sqrt((*per->covariance())(0,0));
1175 if(std::abs(d0/sig_d0) > m_d0_cut) pass = true;
1176 return pass;
1177 }
1178
1179 JpsiXPlusDisplaced::XiCandidate JpsiXPlusDisplaced::getXiCandidate(const xAOD::Vertex* V0vtx, const V0Enum V0, const xAOD::TrackParticle* track3) const {
1180 XiCandidate disVtx;
1181
1182 std::vector<const xAOD::TrackParticle*> tracksV0;
1183 tracksV0.reserve(V0vtx->nTrackParticles());
1184 for(size_t i=0; i<V0vtx->nTrackParticles(); i++) tracksV0.push_back(V0vtx->trackParticle(i));
1185 std::vector<double> massesV0;
1186 if(V0==LAMBDA) massesV0 = m_massesV0_ppi;
1187 else if(V0==LAMBDABAR) massesV0 = m_massesV0_pip;
1188 else if(V0==KS) massesV0 = m_massesV0_pipi;
1189
1190 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
1191 int robustness = 0;
1192 m_iVertexFitter->setRobustness(robustness, *state);
1193 std::vector<Trk::VertexID> vrtList;
1194 // V0 vertex
1195 Trk::VertexID vID = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1196 vrtList.push_back(vID);
1197 // Mother vertex
1198 std::vector<const xAOD::TrackParticle*> tracksDis3{track3};
1199 std::vector<double> massesDis3{m_disVDaug3MassHypo};
1200 m_iVertexFitter->nextVertex(tracksDis3,massesDis3,vrtList,*state);
1201 // Do the work
1202 std::unique_ptr<Trk::VxCascadeInfo> cascade_info = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state) );
1203 if(cascade_info) {
1204 cascade_info->setSVOwnership(true);
1205 double chi2NDF = cascade_info->fitChi2()/cascade_info->nDoF();
1206 if(m_chi2cut_DisV<=0 || chi2NDF < m_chi2cut_DisV) {
1207 const std::vector<std::vector<TLorentzVector> > &moms = cascade_info->getParticleMoms();
1208 const std::vector<xAOD::Vertex*> &cascadeVertices = cascade_info->vertices();
1209 double lxy_SV1_sub = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
1210 if(lxy_SV1_sub > 0.2) {
1211 TLorentzVector totalMom;
1212 for(size_t it=0; it<moms[1].size(); it++) totalMom += moms[1][it];
1213 double disV_mass = totalMom.M();
1214 if(m_useImprovedMass) {
1215 if((V0==LAMBDA || V0==LAMBDABAR) && m_massLd>0) disV_mass += - moms[1][1].M() + m_massLd;
1216 else if(V0==KS && m_massKs>0) disV_mass += - moms[1][1].M() + m_massKs;
1217 }
1218 if(disV_mass>m_DisplacedMassLower && disV_mass<m_DisplacedMassUpper) {
1219 disVtx.track = track3; disVtx.V0vtx = V0vtx; disVtx.V0type = V0;
1220 disVtx.chi2NDF = chi2NDF; disVtx.p4_V0track1 = moms[0][0];
1221 disVtx.p4_V0track2 = moms[0][1]; disVtx.p4_disVtrack = moms[1][0];
1222 }
1223 }
1224 }
1225 }
1226 return disVtx;
1227 }
1228
1229 std::unique_ptr<xAOD::Vertex> JpsiXPlusDisplaced::fitTracks(const xAOD::TrackParticle* track1, const xAOD::TrackParticle* track2, const xAOD::TrackParticle* track3) const {
1230 // Starting point
1231 const Trk::Perigee& aPerigee1 = track1->perigeeParameters();
1232 const Trk::Perigee& aPerigee2 = track2->perigeeParameters();
1233 int sflag(0), errorcode(0);
1234 Amg::Vector3D startingPoint = m_vertexEstimator->getCirclesIntersectionPoint(&aPerigee1,&aPerigee2,sflag,errorcode);
1235 if(errorcode) startingPoint(0) = startingPoint(1) = startingPoint(2) = 0.0;
1236 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
1237 // do the fit
1238 if(track3) {
1239 std::unique_ptr<xAOD::Vertex> fittedVertex( m_iVertexFitter->fit(std::vector<const xAOD::TrackParticle*>{track1,track2,track3}, startingPoint, *state) );
1240 return fittedVertex;
1241 }
1242 else {
1243 std::unique_ptr<xAOD::Vertex> fittedVertex( m_iVertexFitter->fit(std::vector<const xAOD::TrackParticle*>{track1,track2}, startingPoint, *state) );
1244 return fittedVertex;
1245 }
1246 }
1247
1248 JpsiXPlusDisplaced::MesonCandidate JpsiXPlusDisplaced::getDpmCandidate(const xAOD::Vertex* JXvtx, const xAOD::TrackParticle* extraTrk1, const xAOD::TrackParticle* extraTrk2, const xAOD::TrackParticle* extraTrk3) const {
1249 MesonCandidate Dpm;
1250 // Check overlap
1251 std::vector<const xAOD::TrackParticle*> tracksJX;
1252 tracksJX.reserve(JXvtx->nTrackParticles());
1253 for(size_t i=0; i<JXvtx->nTrackParticles(); i++) tracksJX.push_back(JXvtx->trackParticle(i));
1254
1255 TLorentzVector tmp1, tmp2, tmp3;
1256 tmp1.SetPtEtaPhiM(extraTrk1->pt(),extraTrk1->eta(),extraTrk1->phi(),m_extraTrk1MassHypo);
1257 tmp2.SetPtEtaPhiM(extraTrk2->pt(),extraTrk2->eta(),extraTrk2->phi(),m_extraTrk2MassHypo);
1258 tmp3.SetPtEtaPhiM(extraTrk3->pt(),extraTrk3->eta(),extraTrk3->phi(),m_extraTrk3MassHypo);
1259 if((tmp1+tmp2+tmp3).M() < m_DpmMassLower || (tmp1+tmp2+tmp3).M() > m_DpmMassUpper) return Dpm;
1260
1261 std::unique_ptr<xAOD::Vertex> vtx = fitTracks(extraTrk1, extraTrk2, extraTrk3);
1262 if(vtx) {
1263 double chi2NDF = vtx->chiSquared()/vtx->numberDoF();
1264 if(m_chi2cut_Dpm<=0.0 || chi2NDF < m_chi2cut_Dpm) {
1265 double lxyDpm = m_V0Tools->lxy(vtx.get(),JXvtx);
1266 if(lxyDpm>m_lxyDpm_cut) {
1267 Dpm.extraTrack1 = extraTrk1; Dpm.extraTrack2 = extraTrk2; Dpm.extraTrack3 = extraTrk3;
1268 Dpm.chi2NDF = chi2NDF;
1269 TVector3 tot_pt; TVector3 tmp;
1270 for(size_t i=0; i<vtx->vxTrackAtVertex().size(); ++i) {
1271 const Trk::TrackParameters* aPerigee = vtx->vxTrackAtVertex()[i].perigeeAtVertex();
1272 if(aPerigee) {
1273 tmp.SetXYZ(aPerigee->momentum()[Trk::px],aPerigee->momentum()[Trk::py],aPerigee->momentum()[Trk::pz]);
1274 tot_pt += tmp;
1275 }
1276 }
1277 Dpm.pt = tot_pt.Pt();
1278 }
1279 }
1280 }
1281 return Dpm;
1282 }
1283
1284 JpsiXPlusDisplaced::MesonCandidate JpsiXPlusDisplaced::getD0Candidate(const xAOD::Vertex* JXvtx, const xAOD::TrackParticle* extraTrk1, const xAOD::TrackParticle* extraTrk2) const {
1285 MesonCandidate D0;
1286
1287 TLorentzVector tmp1, tmp2;
1288 tmp1.SetPtEtaPhiM(extraTrk1->pt(),extraTrk1->eta(),extraTrk1->phi(),m_extraTrk1MassHypo);
1289 tmp2.SetPtEtaPhiM(extraTrk2->pt(),extraTrk2->eta(),extraTrk2->phi(),m_extraTrk2MassHypo);
1290 if((tmp1+tmp2).M() < m_D0MassLower || (tmp1+tmp2).M() > m_D0MassUpper) return D0;
1291
1292 std::unique_ptr<xAOD::Vertex> vtx = fitTracks(extraTrk1, extraTrk2);
1293 if(vtx) {
1294 double chi2NDF = vtx->chiSquared()/vtx->numberDoF();
1295 if(m_chi2cut_D0<=0.0 || chi2NDF < m_chi2cut_D0) {
1296 double lxyD0 = m_V0Tools->lxy(vtx.get(),JXvtx);
1297 if(lxyD0>m_lxyD0_cut) {
1298 D0.extraTrack1 = extraTrk1; D0.extraTrack2 = extraTrk2;
1299 D0.chi2NDF = chi2NDF;
1300 TVector3 tot_pt; TVector3 tmp;
1301 for(size_t i=0; i<vtx->vxTrackAtVertex().size(); ++i) {
1302 const Trk::TrackParameters* aPerigee = vtx->vxTrackAtVertex()[i].perigeeAtVertex();
1303 if(aPerigee) {
1304 tmp.SetXYZ(aPerigee->momentum()[Trk::px],aPerigee->momentum()[Trk::py],aPerigee->momentum()[Trk::pz]);
1305 tot_pt += tmp;
1306 }
1307 }
1308 D0.pt = tot_pt.Pt();
1309 }
1310 }
1311 }
1312 return D0;
1313 }
1314
1315 std::vector<std::pair<Trk::VxCascadeInfo*,Trk::VxCascadeInfo*> > JpsiXPlusDisplaced::fitMainVtx(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 {
1316 std::vector<std::pair<Trk::VxCascadeInfo*,Trk::VxCascadeInfo*> > result;
1317
1318 std::vector<const xAOD::TrackParticle*> tracksJX;
1319 tracksJX.reserve(JXvtx->nTrackParticles());
1320 for(size_t i=0; i<JXvtx->nTrackParticles(); i++) tracksJX.push_back(JXvtx->trackParticle(i));
1321 if (tracksJX.size() != massesJX.size()) {
1322 ATH_MSG_ERROR("Problems with JX input: number of tracks or track mass inputs is not correct!");
1323 return result;
1324 }
1325 // Check identical tracks in input
1326 if(std::find(tracksJX.cbegin(), tracksJX.cend(), V0vtx->trackParticle(0)) != tracksJX.cend()) return result;
1327 if(std::find(tracksJX.cbegin(), tracksJX.cend(), V0vtx->trackParticle(1)) != tracksJX.cend()) return result;
1328 std::vector<const xAOD::TrackParticle*> tracksV0;
1329 tracksV0.reserve(V0vtx->nTrackParticles());
1330 for(size_t j=0; j<V0vtx->nTrackParticles(); j++) tracksV0.push_back(V0vtx->trackParticle(j));
1331
1332 std::vector<const xAOD::TrackParticle*> tracksJpsi{tracksJX[0], tracksJX[1]};
1333 std::vector<const xAOD::TrackParticle*> tracksX;
1334 if(m_jxDaug_num>=3) tracksX.push_back(tracksJX[2]);
1335 if(m_jxDaug_num==4) tracksX.push_back(tracksJX[3]);
1336
1337 std::vector<double> massesV0;
1338 if(V0==LAMBDA) {
1339 massesV0 = m_massesV0_ppi;
1340 }
1341 else if(V0==LAMBDABAR) {
1342 massesV0 = m_massesV0_pip;
1343 }
1344 else if(V0==KS) {
1345 massesV0 = m_massesV0_pipi;
1346 }
1347
1348 TLorentzVector p4_moth, p4_v0, tmp;
1349 for(size_t it=0; it<JXvtx->nTrackParticles(); it++) {
1350 tmp.SetPtEtaPhiM(JXvtx->trackParticle(it)->pt(), JXvtx->trackParticle(it)->eta(), JXvtx->trackParticle(it)->phi(), massesJX[it]);
1351 p4_moth += tmp;
1352 }
1353 xAOD::BPhysHelper V0_helper(V0vtx);
1354 for(int it=0; it<V0_helper.nRefTrks(); it++) {
1355 p4_moth += V0_helper.refTrk(it,massesV0[it]);
1356 p4_v0 += V0_helper.refTrk(it,massesV0[it]);
1357 }
1358
1359 xAOD::BPhysHelper::pv_type pvtype = xAOD::BPhysHelper::PV_MIN_A0;
1360 if(m_cascadeFitWithPV==1) pvtype = xAOD::BPhysHelper::PV_MAX_SUM_PT2;
1361 else if(m_cascadeFitWithPV==2) pvtype = xAOD::BPhysHelper::PV_MIN_A0;
1362 else if(m_cascadeFitWithPV==3) pvtype = xAOD::BPhysHelper::PV_MIN_Z0;
1363 else if(m_cascadeFitWithPV==4) pvtype = xAOD::BPhysHelper::PV_MIN_Z0_BA;
1364 xAOD::BPhysHelper JX_helper(JXvtx);
1365 const xAOD::Vertex* origPv_xAOD = (m_cascadeFitWithPV>=1 && m_cascadeFitWithPV<=4) ? JX_helper.origPv(pvtype) : nullptr;
1366 const xAOD::Vertex* pv_xAOD = (m_cascadeFitWithPV>=1 && m_cascadeFitWithPV<=4) ? JX_helper.pv(pvtype) : nullptr;
1367 std::unique_ptr<Trk::RecVertex> pv_AOD;
1368 if(pv_xAOD) pv_AOD = std::make_unique<Trk::RecVertex>(pv_xAOD->position(),pv_xAOD->covariancePosition(),pv_xAOD->numberDoF(),pv_xAOD->chiSquared());
1369
1370 SG::AuxElement::Decorator<float> chi2_V1_decor("ChiSquared_V1");
1371 SG::AuxElement::Decorator<int> ndof_V1_decor("nDoF_V1");
1372 SG::AuxElement::Decorator<std::string> type_V1_decor("Type_V1");
1373
1374 SG::AuxElement::Accessor<int> mAcc_gfit("gamma_fit");
1375 SG::AuxElement::Accessor<float> mAcc_gmass("gamma_mass");
1376 SG::AuxElement::Accessor<float> mAcc_gmasserr("gamma_massError");
1377 SG::AuxElement::Accessor<float> mAcc_gchisq("gamma_chisq");
1378 SG::AuxElement::Accessor<int> mAcc_gndof("gamma_ndof");
1379 SG::AuxElement::Accessor<float> mAcc_gprob("gamma_probability");
1380
1381 SG::AuxElement::Decorator<int> mDec_gfit("gamma_fit");
1382 SG::AuxElement::Decorator<float> mDec_gmass("gamma_mass");
1383 SG::AuxElement::Decorator<float> mDec_gmasserr("gamma_massError");
1384 SG::AuxElement::Decorator<float> mDec_gchisq("gamma_chisq");
1385 SG::AuxElement::Decorator<int> mDec_gndof("gamma_ndof");
1386 SG::AuxElement::Decorator<float> mDec_gprob("gamma_probability");
1387 SG::AuxElement::Decorator< std::vector<float> > trk_pxDeco("TrackPx_V0nc");
1388 SG::AuxElement::Decorator< std::vector<float> > trk_pyDeco("TrackPy_V0nc");
1389 SG::AuxElement::Decorator< std::vector<float> > trk_pzDeco("TrackPz_V0nc");
1390
1391 std::vector<float> trk_px;
1392 std::vector<float> trk_py;
1393 std::vector<float> trk_pz;
1394
1395 if(m_extraTrk1MassHypo<=0) {
1396 double main_mass = p4_moth.M();
1397 if(m_useImprovedMass) {
1398 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_v0).M() + m_massJpsi;
1399 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_v0).M() + m_massJX;
1400 if((V0==LAMBDA || V0==LAMBDABAR) && m_massLd>0) main_mass += - p4_v0.M() + m_massLd;
1401 else if(V0==KS && m_massKs>0) main_mass += - p4_v0.M() + m_massKs;
1402 }
1403 if (main_mass < m_MassLower || main_mass > m_MassUpper) return result;
1404
1405 // Apply the user's settings to the fitter
1406 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
1407 // Robustness: http://cdsweb.cern.ch/record/685551
1408 int robustness = 0;
1409 m_iVertexFitter->setRobustness(robustness, *state);
1410 // Build up the topology
1411 // Vertex list
1412 std::vector<Trk::VertexID> vrtList;
1413 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
1414 // V0 vertex
1415 Trk::VertexID vID1;
1416 if (m_constrV0) {
1417 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,V0==KS?m_massKs:m_massLd);
1418 } else {
1419 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1420 }
1421 vrtList.push_back(vID1);
1422 Trk::VertexID vID2;
1423 if(m_JXSubVtx) {
1424 // JX vertex
1425 if (m_constrJX && m_jxDaug_num>2) {
1426 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state,m_massJX);
1427 } else {
1428 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state);
1429 }
1430 vrtList.push_back(vID2);
1431 // Mother vertex including JX and V0
1432 std::vector<const xAOD::TrackParticle*> tp;
1433 std::vector<double> tp_masses;
1434 if(m_constrMainV) {
1435 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList,*state,m_massMainV);
1436 } else {
1437 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList,*state);
1438 }
1439 }
1440 else { // m_JXSubVtx=false
1441 // Mother vertex including JX and V0
1442 if(m_constrMainV) {
1443 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state,m_massMainV);
1444 } else {
1445 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
1446 }
1447 if (m_constrJX && m_jxDaug_num>2) {
1448 std::vector<Trk::VertexID> cnstV;
1449 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
1450 ATH_MSG_WARNING("addMassConstraint for JX failed");
1451 }
1452 }
1453 }
1454 if (m_constrJpsi) {
1455 std::vector<Trk::VertexID> cnstV;
1456 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
1457 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
1458 }
1459 }
1460 if (m_constrX && m_jxDaug_num==4) {
1461 std::vector<Trk::VertexID> cnstV;
1462 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksX,cnstV,*state,m_massX).isSuccess() ) {
1463 ATH_MSG_WARNING("addMassConstraint for X failed");
1464 }
1465 }
1466 // Do the work
1467 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) );
1468
1469 if (fit_result) {
1470 for(auto& v : fit_result->vertices()) {
1471 if(v->nTrackParticles()==0) {
1472 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
1473 v->setTrackParticleLinks(nullLinkVector);
1474 }
1475 }
1476 // reset links to original tracks
1477 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
1478
1479 // necessary to prevent memory leak
1480 fit_result->setSVOwnership(true);
1481
1482 // Chi2/DOF cut
1483 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
1484 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
1485
1486 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
1487 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
1488 size_t iMoth = cascadeVertices.size()-1;
1489 double lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[iMoth]);
1490 if(chi2CutPassed && lxy_SV1>m_lxyV0_cut) {
1491 chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
1492 ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
1493 if(V0==LAMBDA) {
1494 type_V1_decor(*cascadeVertices[0]) = "Lambda";
1495 }
1496 else if(V0==LAMBDABAR) {
1497 type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
1498 }
1499 else if(V0==KS) {
1500 type_V1_decor(*cascadeVertices[0]) = "Ks";
1501 }
1502 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
1503 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
1504 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
1505 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
1506 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
1507 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
1508 trk_px.clear(); trk_py.clear(); trk_pz.clear();
1509 trk_px.reserve(V0_helper.nRefTrks());
1510 trk_py.reserve(V0_helper.nRefTrks());
1511 trk_pz.reserve(V0_helper.nRefTrks());
1512 for(auto&& vec3 : V0_helper.refTrks()) {
1513 trk_px.push_back( vec3.Px() );
1514 trk_py.push_back( vec3.Py() );
1515 trk_pz.push_back( vec3.Pz() );
1516 }
1517 trk_pxDeco(*cascadeVertices[0]) = trk_px;
1518 trk_pyDeco(*cascadeVertices[0]) = trk_py;
1519 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
1520
1521 result.push_back( std::make_pair(fit_result.release(),nullptr) );
1522 }
1523 }
1524 } // for m_extraTrk1MassHypo<=0
1525 else if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo<=0) {
1526 std::vector<double> massesExtra{m_extraTrk1MassHypo};
1527 std::vector<double> massesJXExtra = massesJX; massesJXExtra.push_back(m_extraTrk1MassHypo);
1528
1529 for(const xAOD::TrackParticle* tpExtra : *trackContainer) {
1530 if( tpExtra->pt()<m_extraTrk1MinPt ) continue;
1531 if( !m_trkSelector->decision(*tpExtra, nullptr) ) continue;
1532 // Check identical tracks in input
1533 if(std::find(tracksJX.cbegin(),tracksJX.cend(),tpExtra) != tracksJX.cend()) continue;
1534 if(std::find(tracksV0.cbegin(),tracksV0.cend(),tpExtra) != tracksV0.cend()) continue;
1535
1536 TLorentzVector tmp;
1537 tmp.SetPtEtaPhiM(tpExtra->pt(),tpExtra->eta(),tpExtra->phi(),m_extraTrk1MassHypo);
1538 double main_mass = (p4_moth+tmp).M();
1539 if(m_useImprovedMass) {
1540 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_v0).M() + m_massJpsi;
1541 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_v0).M() + m_massJX;
1542 if((V0==LAMBDA || V0==LAMBDABAR) && m_massLd>0) main_mass += - p4_v0.M() + m_massLd;
1543 else if(V0==KS && m_massKs>0) main_mass += - p4_v0.M() + m_massKs;
1544 }
1545 if(main_mass < m_MassLower || main_mass > m_MassUpper) continue;
1546
1547 std::vector<const xAOD::TrackParticle*> tracksExtra{tpExtra};
1548 std::vector<const xAOD::TrackParticle*> tracksJXExtra = tracksJX; tracksJXExtra.push_back(tpExtra);
1549
1550 // Apply the user's settings to the fitter
1551 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
1552 // Robustness: http://cdsweb.cern.ch/record/685551
1553 int robustness = 0;
1554 m_iVertexFitter->setRobustness(robustness, *state);
1555 // Build up the topology
1556 // Vertex list
1557 std::vector<Trk::VertexID> vrtList;
1558 std::vector<Trk::VertexID> vrtList2;
1559 Trk::VertexID vID2;
1560 if(m_JXV0SubVtx) {
1561 // V0 vertex
1562 Trk::VertexID vID1;
1563 if (m_constrV0) {
1564 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,V0==KS?m_massKs:m_massLd);
1565 } else {
1566 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1567 }
1568 vrtList.push_back(vID1);
1569 // JX+V0 vertex
1570 if(m_constrJXV0) {
1571 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state,m_massJXV0);
1572 } else {
1573 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
1574 }
1575 vrtList2.push_back(vID2);
1576 // Main vertex
1577 if(m_constrMainV) {
1578 m_iVertexFitter->nextVertex(tracksExtra,massesExtra,vrtList2,*state,m_massMainV);
1579 } else {
1580 m_iVertexFitter->nextVertex(tracksExtra,massesExtra,vrtList2,*state);
1581 }
1582 }
1583 else { // m_JXV0SubVtx==false
1584 // V0 vertex
1585 Trk::VertexID vID1;
1586 if (m_constrV0) {
1587 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,V0==KS?m_massKs:m_massLd);
1588 } else {
1589 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1590 }
1591 vrtList.push_back(vID1);
1592 if(m_JXSubVtx) {
1593 // JX vertex
1594 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state);
1595 vrtList.push_back(vID2);
1596 // Mother vertex
1597 if(m_constrMainV) {
1598 m_iVertexFitter->nextVertex(tracksExtra,massesExtra,vrtList,*state,m_massMainV);
1599 } else {
1600 m_iVertexFitter->nextVertex(tracksExtra,massesExtra,vrtList,*state);
1601 }
1602 }
1603 else { // m_JXSubVtx=false
1604 // Mother vertex includes one subvertex (V0) and JX tracks + extra track
1605 if(m_constrMainV) {
1606 vID2 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList,*state,m_massMainV);
1607 } else {
1608 vID2 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList,*state);
1609 }
1610 }
1611 }
1612 if (m_constrJX && m_jxDaug_num>2) {
1613 std::vector<Trk::VertexID> cnstV;
1614 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
1615 ATH_MSG_WARNING("addMassConstraint for JX failed");
1616 }
1617 }
1618 if (m_constrJpsi) {
1619 std::vector<Trk::VertexID> cnstV;
1620 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
1621 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
1622 }
1623 }
1624 if (m_constrX && m_jxDaug_num==4) {
1625 std::vector<Trk::VertexID> cnstV;
1626 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksX,cnstV,*state,m_massX).isSuccess() ) {
1627 ATH_MSG_WARNING("addMassConstraint for X failed");
1628 }
1629 }
1630 // Do the work
1631 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) );
1632
1633 if (fit_result) {
1634 for(auto& v : fit_result->vertices()) {
1635 if(v->nTrackParticles()==0) {
1636 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
1637 v->setTrackParticleLinks(nullLinkVector);
1638 }
1639 }
1640 // reset links to original tracks
1641 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
1642
1643 // necessary to prevent memory leak
1644 fit_result->setSVOwnership(true);
1645
1646 // Chi2/DOF cut
1647 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
1648 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
1649 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
1650 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
1651 size_t iMoth = cascadeVertices.size()-1;
1652 double lxy_SV1(0);
1653 if(m_JXV0SubVtx) {
1654 lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
1655 }
1656 else {
1657 lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[iMoth]);
1658 }
1659 if(chi2CutPassed && lxy_SV1>m_lxyV0_cut) {
1660 chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
1661 ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
1662 if(V0==LAMBDA) {
1663 type_V1_decor(*cascadeVertices[0]) = "Lambda";
1664 }
1665 else if(V0==LAMBDABAR) {
1666 type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
1667 }
1668 else if(V0==KS) {
1669 type_V1_decor(*cascadeVertices[0]) = "Ks";
1670 }
1671 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
1672 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
1673 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
1674 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
1675 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
1676 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
1677 trk_px.clear(); trk_py.clear(); trk_pz.clear();
1678 trk_px.reserve(V0_helper.nRefTrks());
1679 trk_py.reserve(V0_helper.nRefTrks());
1680 trk_pz.reserve(V0_helper.nRefTrks());
1681 for(auto&& vec3 : V0_helper.refTrks()) {
1682 trk_px.push_back( vec3.Px() );
1683 trk_py.push_back( vec3.Py() );
1684 trk_pz.push_back( vec3.Pz() );
1685 }
1686 trk_pxDeco(*cascadeVertices[0]) = trk_px;
1687 trk_pyDeco(*cascadeVertices[0]) = trk_py;
1688 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
1689
1690 result.push_back( std::make_pair(fit_result.release(),nullptr) );
1691 }
1692 }
1693 } // loop over trackContainer
1694 } // for m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo<=0
1695 else if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0 && m_extraTrk3MassHypo<=0) {
1696 std::vector<const xAOD::TrackParticle*> tracksPlus;
1697 std::vector<const xAOD::TrackParticle*> tracksMinus;
1698 for(const xAOD::TrackParticle* tpExtra : *trackContainer) {
1699 if( tpExtra->pt() < std::fmin(m_extraTrk1MinPt,m_extraTrk2MinPt) ) continue;
1700 if( !m_trkSelector->decision(*tpExtra, nullptr) ) continue;
1701 // Check identical tracks in input
1702 if(std::find(tracksJX.cbegin(),tracksJX.cend(),tpExtra) != tracksJX.cend()) continue;
1703 if(std::find(tracksV0.cbegin(),tracksV0.cend(),tpExtra) != tracksV0.cend()) continue;
1704 if(tpExtra->charge()>0) {
1705 tracksPlus.push_back(tpExtra);
1706 }
1707 else {
1708 tracksMinus.push_back(tpExtra);
1709 }
1710 }
1711
1712 MesonCandidateVector D0Candidates(m_maxMesonCandidates, m_MesonPtOrdering);
1713 TLorentzVector p4_ExtraTrk1, p4_ExtraTrk2;
1714 for(const xAOD::TrackParticle* tp1 : tracksPlus) {
1715 for(const xAOD::TrackParticle* tp2 : tracksMinus) {
1716 if((tp1->pt()>m_extraTrk1MinPt && tp2->pt()>m_extraTrk2MinPt) ||
1717 (tp1->pt()>m_extraTrk2MinPt && tp2->pt()>m_extraTrk1MinPt)) {
1718 p4_ExtraTrk1.SetPtEtaPhiM(tp1->pt(), tp1->eta(), tp1->phi(), m_extraTrk1MassHypo);
1719 p4_ExtraTrk2.SetPtEtaPhiM(tp2->pt(), tp2->eta(), tp2->phi(), m_extraTrk2MassHypo);
1720 double main_mass = (p4_moth+p4_ExtraTrk1+p4_ExtraTrk2).M();
1721 if(m_useImprovedMass) {
1722 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_v0).M() + m_massJpsi;
1723 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_v0).M() + m_massJX;
1724 if((V0==LAMBDA || V0==LAMBDABAR) && m_massLd>0) main_mass += - p4_v0.M() + m_massLd;
1725 else if(V0==KS && m_massKs>0) main_mass += - p4_v0.M() + m_massKs;
1726 if(m_massD0>0) main_mass += - (p4_ExtraTrk1+p4_ExtraTrk2).M() + m_massD0;
1727 }
1728 if(main_mass < m_MassLower || main_mass > m_MassUpper) continue;
1729 auto D0 = getD0Candidate(JXvtx,tp1,tp2);
1730 if(D0.extraTrack1) D0Candidates.push_back(D0);
1731 }
1732 }
1733 }
1734
1735 std::vector<double> massesExtra{m_extraTrk1MassHypo,m_extraTrk2MassHypo};
1736
1737 for(auto&& D0 : D0Candidates.vector()) {
1738 std::vector<const xAOD::TrackParticle*> tracksExtra{D0.extraTrack1,D0.extraTrack2};
1739
1740 // Apply the user's settings to the fitter
1741 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
1742 // Robustness: http://cdsweb.cern.ch/record/685551
1743 int robustness = 0;
1744 m_iVertexFitter->setRobustness(robustness, *state);
1745 // Build up the topology
1746 // Vertex list
1747 std::vector<Trk::VertexID> vrtList;
1748 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
1749 // V0 vertex
1750 Trk::VertexID vID1;
1751 if (m_constrV0) {
1752 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,V0==KS?m_massKs:m_massLd);
1753 } else {
1754 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1755 }
1756 vrtList.push_back(vID1);
1757 // D0 vertex
1758 Trk::VertexID vID2;
1759 if (m_constrD0) {
1760 vID2 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state,m_massD0);
1761 } else {
1762 vID2 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state);
1763 }
1764 vrtList.push_back(vID2);
1765 // Mother vertex includes one subvertex (V0) and JX tracks + extra track
1766 Trk::VertexID vID3;
1767 if(m_constrMainV) {
1768 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state,m_massMainV);
1769 } else {
1770 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
1771 }
1772 if (m_constrJX && m_jxDaug_num>2) {
1773 std::vector<Trk::VertexID> cnstV;
1774 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
1775 ATH_MSG_WARNING("addMassConstraint for JX failed");
1776 }
1777 }
1778 if (m_constrJpsi) {
1779 std::vector<Trk::VertexID> cnstV;
1780 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
1781 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
1782 }
1783 }
1784 if (m_constrX && m_jxDaug_num==4) {
1785 std::vector<Trk::VertexID> cnstV;
1786 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksX,cnstV,*state,m_massX).isSuccess() ) {
1787 ATH_MSG_WARNING("addMassConstraint for X failed");
1788 }
1789 }
1790 // Do the work
1791 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) );
1792
1793 if (fit_result) {
1794 for(auto& v : fit_result->vertices()) {
1795 if(v->nTrackParticles()==0) {
1796 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
1797 v->setTrackParticleLinks(nullLinkVector);
1798 }
1799 }
1800 // reset links to original tracks
1801 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
1802
1803 // necessary to prevent memory leak
1804 fit_result->setSVOwnership(true);
1805
1806 // Chi2/DOF cut
1807 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
1808 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
1809 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
1810 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
1811 size_t iMoth = cascadeVertices.size()-1;
1812 double lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[iMoth]);
1813 double lxy_SV2 = m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[iMoth]);
1814 if(chi2CutPassed && lxy_SV1>m_lxyV0_cut && lxy_SV2>m_lxyD0_cut) {
1815 chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
1816 ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
1817 if(V0==LAMBDA) {
1818 type_V1_decor(*cascadeVertices[0]) = "Lambda";
1819 }
1820 else if(V0==LAMBDABAR) {
1821 type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
1822 }
1823 else if(V0==KS) {
1824 type_V1_decor(*cascadeVertices[0]) = "Ks";
1825 }
1826 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
1827 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
1828 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
1829 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
1830 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
1831 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
1832 trk_px.clear(); trk_py.clear(); trk_pz.clear();
1833 trk_px.reserve(V0_helper.nRefTrks());
1834 trk_py.reserve(V0_helper.nRefTrks());
1835 trk_pz.reserve(V0_helper.nRefTrks());
1836 for(auto&& vec3 : V0_helper.refTrks()) {
1837 trk_px.push_back( vec3.Px() );
1838 trk_py.push_back( vec3.Py() );
1839 trk_pz.push_back( vec3.Pz() );
1840 }
1841 trk_pxDeco(*cascadeVertices[0]) = trk_px;
1842 trk_pyDeco(*cascadeVertices[0]) = trk_py;
1843 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
1844
1845 result.push_back( std::make_pair(fit_result.release(),nullptr) );
1846 }
1847 }
1848 } // loop over D0Candidates
1849 } // for m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0 && m_extraTrk3MassHypo<=0
1850 else if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0 && m_extraTrk3MassHypo>0) {
1851 std::vector<const xAOD::TrackParticle*> tracksPlus;
1852 std::vector<const xAOD::TrackParticle*> tracksMinus;
1853 double minTrkPt = std::fmin(std::fmin(m_extraTrk1MinPt,m_extraTrk2MinPt),m_extraTrk3MinPt);
1854 for(const xAOD::TrackParticle* tpExtra : *trackContainer) {
1855 if( tpExtra->pt() < minTrkPt ) continue;
1856 if( !m_trkSelector->decision(*tpExtra, nullptr) ) continue;
1857 // Check identical tracks in input
1858 if(std::find(tracksJX.cbegin(),tracksJX.cend(),tpExtra) != tracksJX.cend()) continue;
1859 if(std::find(tracksV0.cbegin(),tracksV0.cend(),tpExtra) != tracksV0.cend()) continue;
1860 if(tpExtra->charge()>0) {
1861 tracksPlus.push_back(tpExtra);
1862 }
1863 else {
1864 tracksMinus.push_back(tpExtra);
1865 }
1866 }
1867
1868 MesonCandidateVector DpmCandidates(m_maxMesonCandidates, m_MesonPtOrdering);
1869 TLorentzVector p4_ExtraTrk1, p4_ExtraTrk2, p4_ExtraTrk3;
1870 // +,-,- combination (D- -> K+ pi- pi-)
1871 for(const xAOD::TrackParticle* tp1 : tracksPlus) {
1872 if( tp1->pt() < m_extraTrk1MinPt ) continue;
1873 for(auto tp2Itr=tracksMinus.cbegin(); tp2Itr!=tracksMinus.cend(); ++tp2Itr) {
1874 const xAOD::TrackParticle* tp2 = *tp2Itr;
1875 for(auto tp3Itr=tp2Itr+1; tp3Itr!=tracksMinus.cend(); ++tp3Itr) {
1876 const xAOD::TrackParticle* tp3 = *tp3Itr;
1877 if((tp2->pt()>m_extraTrk2MinPt && tp3->pt()>m_extraTrk3MinPt) ||
1878 (tp2->pt()>m_extraTrk3MinPt && tp3->pt()>m_extraTrk2MinPt)) {
1879 p4_ExtraTrk1.SetPtEtaPhiM(tp1->pt(), tp1->eta(), tp1->phi(), m_extraTrk1MassHypo);
1880 p4_ExtraTrk2.SetPtEtaPhiM(tp2->pt(), tp2->eta(), tp2->phi(), m_extraTrk2MassHypo);
1881 p4_ExtraTrk3.SetPtEtaPhiM(tp3->pt(), tp3->eta(), tp3->phi(), m_extraTrk3MassHypo);
1882 double main_mass = (p4_moth+p4_ExtraTrk1+p4_ExtraTrk2+p4_ExtraTrk3).M();
1883 if(m_useImprovedMass) {
1884 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_v0).M() + m_massJpsi;
1885 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_v0).M() + m_massJX;
1886 if((V0==LAMBDA || V0==LAMBDABAR) && m_massLd>0) main_mass += - p4_v0.M() + m_massLd;
1887 else if(V0==KS && m_massKs>0) main_mass += - p4_v0.M() + m_massKs;
1888 if(m_massDpm>0) main_mass += - (p4_ExtraTrk1+p4_ExtraTrk2+p4_ExtraTrk3).M() + m_massDpm;
1889 }
1890 if(main_mass < m_MassLower || main_mass > m_MassUpper) continue;
1891 auto Dpm = getDpmCandidate(JXvtx,tp1,tp2,tp3);
1892 if(Dpm.extraTrack1) DpmCandidates.push_back(Dpm);
1893 }
1894 }
1895 }
1896 }
1897 // -,+,+ combination (D+ -> K- pi+ pi+)
1898 for(const xAOD::TrackParticle* tp1 : tracksMinus) {
1899 if( tp1->pt() < m_extraTrk1MinPt ) continue;
1900 for(auto tp2Itr=tracksPlus.cbegin(); tp2Itr!=tracksPlus.cend(); ++tp2Itr) {
1901 const xAOD::TrackParticle* tp2 = *tp2Itr;
1902 for(auto tp3Itr=tp2Itr+1; tp3Itr!=tracksPlus.cend(); ++tp3Itr) {
1903 const xAOD::TrackParticle* tp3 = *tp3Itr;
1904 if((tp2->pt()>m_extraTrk2MinPt && tp3->pt()>m_extraTrk3MinPt) ||
1905 (tp2->pt()>m_extraTrk3MinPt && tp3->pt()>m_extraTrk2MinPt)) {
1906 p4_ExtraTrk1.SetPtEtaPhiM(tp1->pt(), tp1->eta(), tp1->phi(), m_extraTrk1MassHypo);
1907 p4_ExtraTrk2.SetPtEtaPhiM(tp2->pt(), tp2->eta(), tp2->phi(), m_extraTrk2MassHypo);
1908 p4_ExtraTrk3.SetPtEtaPhiM(tp3->pt(), tp3->eta(), tp3->phi(), m_extraTrk3MassHypo);
1909 double main_mass = (p4_moth+p4_ExtraTrk1+p4_ExtraTrk2+p4_ExtraTrk3).M();
1910 if(m_useImprovedMass) {
1911 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_v0).M() + m_massJpsi;
1912 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_v0).M() + m_massJX;
1913 if((V0==LAMBDA || V0==LAMBDABAR) && m_massLd>0) main_mass += - p4_v0.M() + m_massLd;
1914 else if(V0==KS && m_massKs>0) main_mass += - p4_v0.M() + m_massKs;
1915 if(m_massDpm>0) main_mass += - (p4_ExtraTrk1+p4_ExtraTrk2+p4_ExtraTrk3).M() + m_massDpm;
1916 }
1917 if(main_mass < m_MassLower || main_mass > m_MassUpper) continue;
1918 auto Dpm = getDpmCandidate(JXvtx,tp1,tp2,tp3);
1919 if(Dpm.extraTrack1) DpmCandidates.push_back(Dpm);
1920 }
1921 }
1922 }
1923 }
1924
1925 std::vector<double> massesExtra{m_extraTrk1MassHypo,m_extraTrk2MassHypo,m_extraTrk3MassHypo};
1926
1927 for(auto&& Dpm : DpmCandidates.vector()) {
1928 std::vector<const xAOD::TrackParticle*> tracksExtra{Dpm.extraTrack1,Dpm.extraTrack2,Dpm.extraTrack3};
1929
1930 // Apply the user's settings to the fitter
1931 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
1932 // Robustness: http://cdsweb.cern.ch/record/685551
1933 int robustness = 0;
1934 m_iVertexFitter->setRobustness(robustness, *state);
1935 // Build up the topology
1936 // Vertex list
1937 std::vector<Trk::VertexID> vrtList;
1938 std::vector<Trk::VertexID> vrtList2;
1939 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
1940 if(m_JXV0SubVtx) {
1941 // V0 vertex
1942 Trk::VertexID vID1;
1943 if (m_constrV0) {
1944 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,V0==KS?m_massKs:m_massLd);
1945 } else {
1946 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1947 }
1948 vrtList.push_back(vID1);
1949 // JX+V0 vertex
1950 Trk::VertexID vID2;
1951 if(m_constrJXV0) {
1952 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state,m_massJXV0);
1953 } else {
1954 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
1955 }
1956 vrtList2.push_back(vID2);
1957 Trk::VertexID vID3;
1958 if (m_constrDpm) {
1959 vID3 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state,m_massDpm);
1960 } else {
1961 vID3 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state);
1962 }
1963 vrtList2.push_back(vID3);
1964 // Mother vertex includes two subvertices: V0, Dpm and JX tracks
1965 std::vector<const xAOD::TrackParticle*> tp;
1966 std::vector<double> tp_masses;
1967 if(m_constrMainV) {
1968 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2,*state,m_massMainV);
1969 } else {
1970 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2,*state);
1971 }
1972 if (m_constrJX && m_jxDaug_num>2) {
1973 std::vector<Trk::VertexID> cnstV;
1974 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
1975 ATH_MSG_WARNING("addMassConstraint for JX failed");
1976 }
1977 }
1978 if (m_constrJpsi) {
1979 std::vector<Trk::VertexID> cnstV;
1980 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
1981 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
1982 }
1983 }
1984 if (m_constrX && m_jxDaug_num==4) {
1985 std::vector<Trk::VertexID> cnstV;
1986 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksX,cnstV,*state,m_massX).isSuccess() ) {
1987 ATH_MSG_WARNING("addMassConstraint for X failed");
1988 }
1989 }
1990 }
1991 else { // m_JXV0SubVtx==false
1992 // V0 vertex
1993 Trk::VertexID vID1;
1994 if (m_constrV0) {
1995 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,V0==KS?m_massKs:m_massLd);
1996 } else {
1997 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1998 }
1999 vrtList.push_back(vID1);
2000 // Dpm vertex
2001 Trk::VertexID vID2;
2002 if (m_constrDpm) {
2003 vID2 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state,m_massDpm);
2004 } else {
2005 vID2 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state);
2006 }
2007 vrtList.push_back(vID2);
2008 // Mother vertex includes two subvertices: V0, Dpm and JX tracks
2009 Trk::VertexID vID3;
2010 if(m_constrMainV) {
2011 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state,m_massMainV);
2012 } else {
2013 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
2014 }
2015 if (m_constrJX && m_jxDaug_num>2) {
2016 std::vector<Trk::VertexID> cnstV;
2017 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
2018 ATH_MSG_WARNING("addMassConstraint for JX failed");
2019 }
2020 }
2021 if (m_constrJpsi) {
2022 std::vector<Trk::VertexID> cnstV;
2023 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
2024 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
2025 }
2026 }
2027 if (m_constrX && m_jxDaug_num==4) {
2028 std::vector<Trk::VertexID> cnstV;
2029 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksX,cnstV,*state,m_massX).isSuccess() ) {
2030 ATH_MSG_WARNING("addMassConstraint for X failed");
2031 }
2032 }
2033 }
2034 // Do the work
2035 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) );
2036
2037 if (fit_result) {
2038 for(auto& v : fit_result->vertices()) {
2039 if(v->nTrackParticles()==0) {
2040 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
2041 v->setTrackParticleLinks(nullLinkVector);
2042 }
2043 }
2044 // reset links to original tracks
2045 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
2046
2047 // necessary to prevent memory leak
2048 fit_result->setSVOwnership(true);
2049
2050 // Chi2/DOF cut
2051 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
2052 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
2053 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
2054 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
2055 size_t iMoth = cascadeVertices.size()-1;
2056 double lxy_SV1(0), lxy_SV2(0);
2057 if(m_JXV0SubVtx) {
2058 lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
2059 lxy_SV2 = m_CascadeTools->lxy(moms[2],cascadeVertices[2],cascadeVertices[iMoth]);
2060 }
2061 else {
2062 lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[iMoth]);
2063 lxy_SV2 = m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[iMoth]);
2064 }
2065 if(chi2CutPassed && lxy_SV1>m_lxyV0_cut && lxy_SV2>m_lxyDpm_cut) {
2066 chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
2067 ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
2068 if(V0==LAMBDA) {
2069 type_V1_decor(*cascadeVertices[0]) = "Lambda";
2070 }
2071 else if(V0==LAMBDABAR) {
2072 type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
2073 }
2074 else if(V0==KS) {
2075 type_V1_decor(*cascadeVertices[0]) = "Ks";
2076 }
2077 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
2078 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
2079 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
2080 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
2081 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
2082 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
2083 trk_px.clear(); trk_py.clear(); trk_pz.clear();
2084 trk_px.reserve(V0_helper.nRefTrks());
2085 trk_py.reserve(V0_helper.nRefTrks());
2086 trk_pz.reserve(V0_helper.nRefTrks());
2087 for(auto&& vec3 : V0_helper.refTrks()) {
2088 trk_px.push_back( vec3.Px() );
2089 trk_py.push_back( vec3.Py() );
2090 trk_pz.push_back( vec3.Pz() );
2091 }
2092 trk_pxDeco(*cascadeVertices[0]) = trk_px;
2093 trk_pyDeco(*cascadeVertices[0]) = trk_py;
2094 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
2095
2096 // refit with main vertex mass constraint if required
2097 if(!m_constrMainV && m_doPostMainVContrFit) {
2098 TLorentzVector totalMom;
2099 for(size_t it=0; it<moms[iMoth].size(); it++) totalMom += moms[iMoth][it];
2100 double mainV_mass = totalMom.M();
2101 if(mainV_mass>m_PostMassLower && mainV_mass<m_PostMassUpper) {
2102 std::unique_ptr<Trk::IVKalState> state_mvc = m_iVertexFitter->makeState();
2103 int robustness_mvc = 0;
2104 m_iVertexFitter->setRobustness(robustness_mvc, *state_mvc);
2105 std::vector<Trk::VertexID> vrtList_mvc;
2106 std::vector<Trk::VertexID> vrtList2_mvc;
2107 if(m_JXV0SubVtx) {
2108 Trk::VertexID vID1_mvc; // V0 vertex
2109 if (m_constrV0) {
2110 vID1_mvc = m_iVertexFitter->startVertex(tracksV0,massesV0,*state_mvc,V0==KS?m_massKs:m_massLd);
2111 } else {
2112 vID1_mvc = m_iVertexFitter->startVertex(tracksV0,massesV0,*state_mvc);
2113 }
2114 vrtList_mvc.push_back(vID1_mvc);
2115 Trk::VertexID vID2_mvc; // JX+V0 vertex
2116 if(m_constrJXV0) {
2117 vID2_mvc = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList_mvc,*state_mvc,m_massJXV0);
2118 } else {
2119 vID2_mvc = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList_mvc,*state_mvc);
2120 }
2121 vrtList2_mvc.push_back(vID2_mvc);
2122 Trk::VertexID vID3_mvc; // Dpm vertex
2123 if (m_constrDpm) {
2124 vID3_mvc = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state_mvc,m_massDpm);
2125 } else {
2126 vID3_mvc = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state_mvc);
2127 }
2128 vrtList2_mvc.push_back(vID3_mvc);
2129 std::vector<const xAOD::TrackParticle*> tp;
2130 std::vector<double> tp_masses;
2131 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2_mvc,*state_mvc,m_massMainV); // Mother vertex
2132 if (m_constrJX && m_jxDaug_num>2) {
2133 std::vector<Trk::VertexID> cnstV_mvc;
2134 if ( !m_iVertexFitter->addMassConstraint(vID2_mvc,tracksJX,cnstV_mvc,*state_mvc,m_massJX).isSuccess() ) {
2135 ATH_MSG_WARNING("addMassConstraint for JX failed");
2136 }
2137 }
2138 if (m_constrJpsi) {
2139 std::vector<Trk::VertexID> cnstV_mvc;
2140 if ( !m_iVertexFitter->addMassConstraint(vID2_mvc,tracksJpsi,cnstV_mvc,*state_mvc,m_massJpsi).isSuccess() ) {
2141 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
2142 }
2143 }
2144 if (m_constrX && m_jxDaug_num==4) {
2145 std::vector<Trk::VertexID> cnstV_mvc;
2146 if ( !m_iVertexFitter->addMassConstraint(vID2_mvc,tracksX,cnstV_mvc,*state_mvc,m_massX).isSuccess() ) {
2147 ATH_MSG_WARNING("addMassConstraint for X failed");
2148 }
2149 }
2150 }
2151 else {
2152 Trk::VertexID vID1_mvc; // V0 vertex
2153 if (m_constrV0) {
2154 vID1_mvc = m_iVertexFitter->startVertex(tracksV0,massesV0,*state_mvc,V0==KS?m_massKs:m_massLd);
2155 } else {
2156 vID1_mvc = m_iVertexFitter->startVertex(tracksV0,massesV0,*state_mvc);
2157 }
2158 vrtList_mvc.push_back(vID1_mvc);
2159 Trk::VertexID vID2_mvc; // Dpm vertex
2160 if (m_constrDpm) {
2161 vID2_mvc = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state_mvc,m_massDpm);
2162 } else {
2163 vID2_mvc = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state_mvc);
2164 }
2165 vrtList_mvc.push_back(vID2_mvc);
2166 Trk::VertexID vID3_mvc = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList_mvc,*state_mvc,m_massMainV); // Mother vertex
2167 if (m_constrJX && m_jxDaug_num>2) {
2168 std::vector<Trk::VertexID> cnstV_mvc;
2169 if ( !m_iVertexFitter->addMassConstraint(vID3_mvc,tracksJX,cnstV_mvc,*state_mvc,m_massJX).isSuccess() ) {
2170 ATH_MSG_WARNING("addMassConstraint for JX failed");
2171 }
2172 }
2173 if (m_constrJpsi) {
2174 std::vector<Trk::VertexID> cnstV_mvc;
2175 if ( !m_iVertexFitter->addMassConstraint(vID3_mvc,tracksJpsi,cnstV_mvc,*state_mvc,m_massJpsi).isSuccess() ) {
2176 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
2177 }
2178 }
2179 if (m_constrX && m_jxDaug_num==4) {
2180 std::vector<Trk::VertexID> cnstV_mvc;
2181 if ( !m_iVertexFitter->addMassConstraint(vID3_mvc,tracksX,cnstV_mvc,*state_mvc,m_massX).isSuccess() ) {
2182 ATH_MSG_WARNING("addMassConstraint for X failed");
2183 }
2184 }
2185 }
2186 // Do the work
2187 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) );
2188
2189 if (fit_result_mvc) {
2190 for(auto& v : fit_result_mvc->vertices()) {
2191 if(v->nTrackParticles()==0) {
2192 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
2193 v->setTrackParticleLinks(nullLinkVector);
2194 }
2195 }
2196 BPhysPVCascadeTools::PrepareVertexLinks(fit_result_mvc.get(), trackCols);
2197 fit_result_mvc->setSVOwnership(true);
2198 chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
2199 ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
2200 if(V0==LAMBDA) {
2201 type_V1_decor(*cascadeVertices[0]) = "Lambda";
2202 }
2203 else if(V0==LAMBDABAR) {
2204 type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
2205 }
2206 else if(V0==KS) {
2207 type_V1_decor(*cascadeVertices[0]) = "Ks";
2208 }
2209 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
2210 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
2211 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
2212 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
2213 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
2214 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
2215 trk_px.clear(); trk_py.clear(); trk_pz.clear();
2216 trk_px.reserve(V0_helper.nRefTrks());
2217 trk_py.reserve(V0_helper.nRefTrks());
2218 trk_pz.reserve(V0_helper.nRefTrks());
2219 for(auto&& vec3 : V0_helper.refTrks()) {
2220 trk_px.push_back( vec3.Px() );
2221 trk_py.push_back( vec3.Py() );
2222 trk_pz.push_back( vec3.Pz() );
2223 }
2224 trk_pxDeco(*cascadeVertices[0]) = trk_px;
2225 trk_pyDeco(*cascadeVertices[0]) = trk_py;
2226 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
2227
2228 result.push_back( std::make_pair(fit_result.release(),fit_result_mvc.release()) );
2229 }
2230 else result.push_back( std::make_pair(fit_result.release(),nullptr) );
2231 }
2232 else result.push_back( std::make_pair(fit_result.release(),nullptr) );
2233 }
2234 else result.push_back( std::make_pair(fit_result.release(),nullptr) );
2235 }
2236 }
2237 } // loop over DpmCandidates
2238 } // for m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0 && m_extraTrk3MassHypo>0
2239
2240 if(pv_xAOD) {
2241 for(auto cascade_info_pair : result) {
2242 if(cascade_info_pair.first && cascade_info_pair.first->getParticleMoms().size()>0) {
2243 size_t index = cascade_info_pair.first->getParticleMoms().size() - 1;
2244 const std::vector<TLorentzVector> &mom = cascade_info_pair.first->getParticleMoms()[index];
2245 const Amg::MatrixX &cov = cascade_info_pair.first->getCovariance()[index];
2246 const xAOD::Vertex* mainVertex = cascade_info_pair.first->vertices()[index];
2247 xAOD::BPhysHypoHelper vtx(m_hypoName, mainVertex);
2248 bool isInDefaultPVCont = false;
2249 for(const xAOD::Vertex* pvVtx : *defaultPVContainer) {
2250 if(pv_xAOD == pvVtx) { isInDefaultPVCont = true; break; }
2251 }
2252 if(isInDefaultPVCont) vtx.setPv( pv_xAOD, defaultPVContainer, pvtype );
2253 else vtx.setPv( pv_xAOD, pvContainer, pvtype );
2254 if(origPv_xAOD) vtx.setOrigPv( origPv_xAOD, defaultPVContainer, pvtype );
2255 vtx.setLxy ( m_CascadeTools->lxy (mom, vtx.vtx(), pv_xAOD), pvtype );
2256 vtx.setLxyErr ( m_CascadeTools->lxyError (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2257 vtx.setA0 ( m_CascadeTools->a0 (mom, vtx.vtx(), pv_xAOD), pvtype );
2258 vtx.setA0Err ( m_CascadeTools->a0Error (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2259 vtx.setA0xy ( m_CascadeTools->a0xy (mom, vtx.vtx(), pv_xAOD), pvtype );
2260 vtx.setA0xyErr( m_CascadeTools->a0xyError(mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2261 vtx.setZ0 ( m_CascadeTools->a0z (mom, vtx.vtx(), pv_xAOD), pvtype );
2262 vtx.setZ0Err ( m_CascadeTools->a0zError (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2263 vtx.setRefitPVStatus( 0, pvtype );
2264 // Proper decay times
2265 vtx.setTau( m_CascadeTools->tau(mom, vtx.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
2266 vtx.setTauErr( m_CascadeTools->tauError(mom, cov, vtx.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
2267 vtx.setTau( m_CascadeTools->tau(mom, vtx.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
2268 vtx.setTauErr( m_CascadeTools->tauError(mom, cov, vtx.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
2269
2270 if(cascade_info_pair.second && cascade_info_pair.second->getParticleMoms().size()>0) {
2271 index = cascade_info_pair.second->getParticleMoms().size() - 1;
2272 const std::vector<TLorentzVector> &mom_mvc = cascade_info_pair.second->getParticleMoms()[index];
2273 const Amg::MatrixX &cov_mvc = cascade_info_pair.second->getCovariance()[index];
2274 const xAOD::Vertex* mainVertex_mvc = cascade_info_pair.second->vertices()[index];
2275 xAOD::BPhysHypoHelper vtx_mvc(m_hypoName, mainVertex_mvc);
2276 if(isInDefaultPVCont) vtx_mvc.setPv( pv_xAOD, defaultPVContainer, pvtype );
2277 else vtx_mvc.setPv( pv_xAOD, pvContainer, pvtype );
2278 if(origPv_xAOD) vtx.setOrigPv( origPv_xAOD, defaultPVContainer, pvtype );
2279 vtx_mvc.setLxy ( m_CascadeTools->lxy (mom_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2280 vtx_mvc.setLxyErr ( m_CascadeTools->lxyError (mom_mvc, cov_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2281 vtx_mvc.setA0 ( m_CascadeTools->a0 (mom_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2282 vtx_mvc.setA0Err ( m_CascadeTools->a0Error (mom_mvc, cov_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2283 vtx_mvc.setA0xy ( m_CascadeTools->a0xy (mom_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2284 vtx_mvc.setA0xyErr( m_CascadeTools->a0xyError(mom_mvc, cov_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2285 vtx_mvc.setZ0 ( m_CascadeTools->a0z (mom_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2286 vtx_mvc.setZ0Err ( m_CascadeTools->a0zError (mom_mvc, cov_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2287 vtx_mvc.setRefitPVStatus( 0, pvtype );
2288 // Proper decay times
2289 vtx_mvc.setTau( m_CascadeTools->tau(mom_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
2290 vtx_mvc.setTauErr( m_CascadeTools->tauError(mom_mvc, cov_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
2291 vtx_mvc.setTau( m_CascadeTools->tau(mom_mvc, vtx_mvc.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
2292 vtx_mvc.setTauErr( m_CascadeTools->tauError(mom_mvc, cov_mvc, vtx_mvc.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
2293 }
2294 }
2295 }
2296 }
2297
2298 return result;
2299 }
2300
2301 std::vector<std::pair<Trk::VxCascadeInfo*,Trk::VxCascadeInfo*> > JpsiXPlusDisplaced::fitMainVtx(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 {
2302 std::vector<std::pair<Trk::VxCascadeInfo*,Trk::VxCascadeInfo*> > result;
2303
2304 std::vector<const xAOD::TrackParticle*> tracksJX;
2305 tracksJX.reserve(JXvtx->nTrackParticles());
2306 for(size_t i=0; i<JXvtx->nTrackParticles(); i++) tracksJX.push_back(JXvtx->trackParticle(i));
2307 if (tracksJX.size() != massesJX.size()) {
2308 ATH_MSG_ERROR("Problems with JX input: number of tracks or track mass inputs is not correct!");
2309 return result;
2310 }
2311 // Check identical tracks in input
2312 if(std::find(tracksJX.cbegin(), tracksJX.cend(), disVtx.V0vtx->trackParticle(0)) != tracksJX.cend()) return result;
2313 if(std::find(tracksJX.cbegin(), tracksJX.cend(), disVtx.V0vtx->trackParticle(1)) != tracksJX.cend()) return result;
2314 std::vector<const xAOD::TrackParticle*> tracksV0;
2315 tracksV0.reserve(disVtx.V0vtx->nTrackParticles());
2316 for(size_t j=0; j<disVtx.V0vtx->nTrackParticles(); j++) tracksV0.push_back(disVtx.V0vtx->trackParticle(j));
2317
2318 if(std::find(tracksJX.cbegin(), tracksJX.cend(), disVtx.track) != tracksJX.cend()) return result;
2319 std::vector<const xAOD::TrackParticle*> tracks3{disVtx.track};
2320 std::vector<double> massesDis3{m_disVDaug3MassHypo};
2321
2322 std::vector<const xAOD::TrackParticle*> tracksJpsi{tracksJX[0], tracksJX[1]};
2323 std::vector<const xAOD::TrackParticle*> tracksX;
2324 if(m_jxDaug_num>=3) tracksX.push_back(tracksJX[2]);
2325 if(m_jxDaug_num==4) tracksX.push_back(tracksJX[3]);
2326
2327 std::vector<double> massesV0;
2328 if(disVtx.V0type==LAMBDA) {
2329 massesV0 = m_massesV0_ppi;
2330 }
2331 else if(disVtx.V0type==LAMBDABAR) {
2332 massesV0 = m_massesV0_pip;
2333 }
2334 else if(disVtx.V0type==KS) {
2335 massesV0 = m_massesV0_pipi;
2336 }
2337
2338 std::vector<double> massesDisV = massesV0; massesDisV.push_back(m_disVDaug3MassHypo);
2339
2340 TLorentzVector p4_moth, p4_disV, tmp;
2341 for(size_t it=0; it<JXvtx->nTrackParticles(); it++) {
2342 tmp.SetPtEtaPhiM(JXvtx->trackParticle(it)->pt(), JXvtx->trackParticle(it)->eta(), JXvtx->trackParticle(it)->phi(), massesJX[it]);
2343 p4_moth += tmp;
2344 }
2345 p4_disV += disVtx.p4_V0track1; p4_disV += disVtx.p4_V0track2; p4_disV += disVtx.p4_disVtrack;
2346 p4_moth += p4_disV;
2347
2348 xAOD::BPhysHelper::pv_type pvtype = xAOD::BPhysHelper::PV_MIN_A0;
2349 if(m_cascadeFitWithPV==1) pvtype = xAOD::BPhysHelper::PV_MAX_SUM_PT2;
2350 else if(m_cascadeFitWithPV==2) pvtype = xAOD::BPhysHelper::PV_MIN_A0;
2351 else if(m_cascadeFitWithPV==3) pvtype = xAOD::BPhysHelper::PV_MIN_Z0;
2352 else if(m_cascadeFitWithPV==4) pvtype = xAOD::BPhysHelper::PV_MIN_Z0_BA;
2353 xAOD::BPhysHelper JX_helper(JXvtx);
2354 const xAOD::Vertex* origPv_xAOD = (m_cascadeFitWithPV>=1 && m_cascadeFitWithPV<=4) ? JX_helper.origPv(pvtype) : nullptr;
2355 const xAOD::Vertex* pv_xAOD = (m_cascadeFitWithPV>=1 && m_cascadeFitWithPV<=4) ? JX_helper.pv(pvtype) : nullptr;
2356 std::unique_ptr<Trk::RecVertex> pv_AOD;
2357 if(pv_xAOD) pv_AOD = std::make_unique<Trk::RecVertex>(pv_xAOD->position(),pv_xAOD->covariancePosition(),pv_xAOD->numberDoF(),pv_xAOD->chiSquared());
2358
2359 SG::AuxElement::Decorator<float> chi2_V0_decor("ChiSquared_V0");
2360 SG::AuxElement::Decorator<int> ndof_V0_decor("nDoF_V0");
2361 SG::AuxElement::Decorator<std::string> type_V0_decor("Type_V0");
2362
2363 SG::AuxElement::Accessor<int> mAcc_gfit("gamma_fit");
2364 SG::AuxElement::Accessor<float> mAcc_gmass("gamma_mass");
2365 SG::AuxElement::Accessor<float> mAcc_gmasserr("gamma_massError");
2366 SG::AuxElement::Accessor<float> mAcc_gchisq("gamma_chisq");
2367 SG::AuxElement::Accessor<int> mAcc_gndof("gamma_ndof");
2368 SG::AuxElement::Accessor<float> mAcc_gprob("gamma_probability");
2369
2370 SG::AuxElement::Decorator<int> mDec_gfit("gamma_fit");
2371 SG::AuxElement::Decorator<float> mDec_gmass("gamma_mass");
2372 SG::AuxElement::Decorator<float> mDec_gmasserr("gamma_massError");
2373 SG::AuxElement::Decorator<float> mDec_gchisq("gamma_chisq");
2374 SG::AuxElement::Decorator<int> mDec_gndof("gamma_ndof");
2375 SG::AuxElement::Decorator<float> mDec_gprob("gamma_probability");
2376 SG::AuxElement::Decorator< std::vector<float> > trk_pxDeco("TrackPx_V0nc");
2377 SG::AuxElement::Decorator< std::vector<float> > trk_pyDeco("TrackPy_V0nc");
2378 SG::AuxElement::Decorator< std::vector<float> > trk_pzDeco("TrackPz_V0nc");
2379 SG::AuxElement::Decorator<float> trk_px_deco("TrackPx_DisVnc");
2380 SG::AuxElement::Decorator<float> trk_py_deco("TrackPy_DisVnc");
2381 SG::AuxElement::Decorator<float> trk_pz_deco("TrackPz_DisVnc");
2382
2383 std::vector<float> trk_px;
2384 std::vector<float> trk_py;
2385 std::vector<float> trk_pz;
2386
2387 if(m_extraTrk1MassHypo<=0) {
2388 double main_mass = p4_moth.M();
2389 if(m_useImprovedMass) {
2390 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_disV).M() + m_massJpsi;
2391 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_disV).M() + m_massJX;
2392 if(m_massDisV>0) main_mass += - p4_disV.M() + m_massDisV;
2393 }
2394 if (main_mass < m_MassLower || main_mass > m_MassUpper) return result;
2395
2396 // Apply the user's settings to the fitter
2397 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
2398 // Robustness: http://cdsweb.cern.ch/record/685551
2399 int robustness = 0;
2400 m_iVertexFitter->setRobustness(robustness, *state);
2401 // Build up the topology
2402 // Vertex list
2403 std::vector<Trk::VertexID> vrtList;
2404 std::vector<Trk::VertexID> vrtList2;
2405 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
2406 // V0 vertex
2407 Trk::VertexID vID1;
2408 if (m_constrV0) {
2409 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,disVtx.V0type==KS?m_massKs:m_massLd);
2410 } else {
2411 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
2412 }
2413 vrtList.push_back(vID1);
2414 // Displaced vertex
2415 Trk::VertexID vID2;
2416 if (m_constrDisV) {
2417 vID2 = m_iVertexFitter->nextVertex(tracks3,massesDis3,vrtList,*state,m_massDisV);
2418 } else {
2419 vID2 = m_iVertexFitter->nextVertex(tracks3,massesDis3,vrtList,*state);
2420 }
2421 vrtList2.push_back(vID2);
2422 Trk::VertexID vID3;
2423 if(m_JXSubVtx) {
2424 // JX vertex
2425 if (m_constrJX && m_jxDaug_num>2) {
2426 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state,m_massJX);
2427 } else {
2428 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state);
2429 }
2430 vrtList2.push_back(vID3);
2431 // Mother vertex includes two subvertices: DisV and JX
2432 std::vector<const xAOD::TrackParticle*> tp;
2433 std::vector<double> tp_masses;
2434 if(m_constrMainV) {
2435 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2,*state,m_massMainV);
2436 } else {
2437 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2,*state);
2438 }
2439 }
2440 else { // m_JXSubVtx=false
2441 // Mother vertex includes just one subvertex (DisV) and JX tracks
2442 if(m_constrMainV) {
2443 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList2,*state,m_massMainV);
2444 } else {
2445 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList2,*state);
2446 }
2447 if (m_constrJX && m_jxDaug_num>2) {
2448 std::vector<Trk::VertexID> cnstV;
2449 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
2450 ATH_MSG_WARNING("addMassConstraint for JX failed");
2451 }
2452 }
2453 }
2454 if (m_constrJpsi) {
2455 std::vector<Trk::VertexID> cnstV;
2456 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
2457 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
2458 }
2459 }
2460 if (m_constrX && m_jxDaug_num==4) {
2461 std::vector<Trk::VertexID> cnstV;
2462 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksX,cnstV,*state,m_massX).isSuccess() ) {
2463 ATH_MSG_WARNING("addMassConstraint for X failed");
2464 }
2465 }
2466 // Do the work
2467 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) );
2468
2469 if (fit_result) {
2470 for(auto& v : fit_result->vertices()) {
2471 if(v->nTrackParticles()==0) {
2472 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
2473 v->setTrackParticleLinks(nullLinkVector);
2474 }
2475 }
2476 // reset links to original tracks
2477 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
2478
2479 // necessary to prevent memory leak
2480 fit_result->setSVOwnership(true);
2481
2482 // Chi2/DOF cut
2483 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
2484 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
2485
2486 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
2487 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
2488 size_t iMoth = cascadeVertices.size()-1;
2489 double lxy_SV1_sub = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
2490 double lxy_SV1 = m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[iMoth]);
2491
2492 if(chi2CutPassed && lxy_SV1>m_lxyDisV_cut && lxy_SV1_sub>m_lxyV0_cut) {
2493 chi2_V0_decor(*cascadeVertices[0]) = disVtx.V0vtx->chiSquared();
2494 ndof_V0_decor(*cascadeVertices[0]) = disVtx.V0vtx->numberDoF();
2495 if(disVtx.V0type==LAMBDA) {
2496 type_V0_decor(*cascadeVertices[0]) = "Lambda";
2497 }
2498 else if(disVtx.V0type==LAMBDABAR) {
2499 type_V0_decor(*cascadeVertices[0]) = "Lambdabar";
2500 }
2501 else if(disVtx.V0type==KS) {
2502 type_V0_decor(*cascadeVertices[0]) = "Ks";
2503 }
2504 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*disVtx.V0vtx) ? mAcc_gfit(*disVtx.V0vtx) : 0;
2505 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*disVtx.V0vtx) ? mAcc_gmass(*disVtx.V0vtx) : -1;
2506 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*disVtx.V0vtx) ? mAcc_gmasserr(*disVtx.V0vtx) : -1;
2507 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*disVtx.V0vtx) ? mAcc_gchisq(*disVtx.V0vtx) : 999999;
2508 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*disVtx.V0vtx) ? mAcc_gndof(*disVtx.V0vtx) : 0;
2509 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*disVtx.V0vtx) ? mAcc_gprob(*disVtx.V0vtx) : -1;
2510 trk_px.clear(); trk_py.clear(); trk_pz.clear();
2511 trk_px.push_back( disVtx.p4_V0track1.Px() ); trk_px.push_back( disVtx.p4_V0track2.Px() );
2512 trk_py.push_back( disVtx.p4_V0track1.Py() ); trk_py.push_back( disVtx.p4_V0track2.Py() );
2513 trk_pz.push_back( disVtx.p4_V0track1.Pz() ); trk_pz.push_back( disVtx.p4_V0track2.Pz() );
2514 trk_pxDeco(*cascadeVertices[0]) = trk_px;
2515 trk_pyDeco(*cascadeVertices[0]) = trk_py;
2516 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
2517 trk_px_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Px();
2518 trk_py_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Py();
2519 trk_pz_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Pz();
2520
2521 result.push_back( std::make_pair(fit_result.release(),nullptr) );
2522 }
2523 }
2524 } // m_extraTrk1MassHypo<=0
2525 else { // m_extraTrk1MassHypo>0
2526 std::vector<double> massesExtra{m_extraTrk1MassHypo};
2527 std::vector<double> massesJXExtra = massesJX; massesJXExtra.push_back(m_extraTrk1MassHypo);
2528
2529 for(const xAOD::TrackParticle* tpExtra : *trackContainer) {
2530 if ( tpExtra->pt()<m_extraTrk1MinPt ) continue;
2531 if ( !m_trkSelector->decision(*tpExtra, nullptr) ) continue;
2532 // Check identical tracks in input
2533 if(std::find(tracksJX.cbegin(),tracksJX.cend(),tpExtra) != tracksJX.cend()) continue;
2534 if(std::find(tracksV0.cbegin(),tracksV0.cend(),tpExtra) != tracksV0.cend()) continue;
2535 if(tpExtra == disVtx.track) continue;
2536
2537 TLorentzVector tmp;
2538 tmp.SetPtEtaPhiM(tpExtra->pt(),tpExtra->eta(),tpExtra->phi(),m_extraTrk1MassHypo);
2539 double main_mass = (p4_moth+tmp).M();
2540 if(m_useImprovedMass) {
2541 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_disV).M() + m_massJpsi;
2542 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_disV).M() + m_massJX;
2543 if(m_massDisV>0) main_mass += - p4_disV.M() + m_massDisV;
2544 }
2545 if (main_mass < m_MassLower || main_mass > m_MassUpper) continue;
2546
2547 std::vector<const xAOD::TrackParticle*> tracksExtra{tpExtra};
2548 std::vector<const xAOD::TrackParticle*> tracksJXExtra = tracksJX; tracksJXExtra.push_back(tpExtra);
2549
2550 // Apply the user's settings to the fitter
2551 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
2552 // Robustness: http://cdsweb.cern.ch/record/685551
2553 int robustness = 0;
2554 m_iVertexFitter->setRobustness(robustness, *state);
2555 // Build up the topology
2556 // Vertex list
2557 std::vector<Trk::VertexID> vrtList;
2558 std::vector<Trk::VertexID> vrtList2;
2559 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
2560 // V0 vertex
2561 Trk::VertexID vID1;
2562 if (m_constrV0) {
2563 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,disVtx.V0type==KS?m_massKs:m_massLd);
2564 } else {
2565 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
2566 }
2567 vrtList.push_back(vID1);
2568 // Displaced vertex
2569 Trk::VertexID vID2;
2570 if (m_constrDisV) {
2571 vID2 = m_iVertexFitter->nextVertex(tracks3,massesDis3,vrtList,*state,m_massDisV);
2572 } else {
2573 vID2 = m_iVertexFitter->nextVertex(tracks3,massesDis3,vrtList,*state);
2574 }
2575 vrtList2.push_back(vID2);
2576 Trk::VertexID vID3;
2577 if(m_JXSubVtx) {
2578 // JXExtra vertex
2579 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state);
2580 vrtList2.push_back(vID3);
2581 // Mother vertex
2582 if(m_constrMainV) {
2583 m_iVertexFitter->nextVertex(tracksExtra,massesExtra,vrtList2,*state,m_massMainV);
2584 } else {
2585 m_iVertexFitter->nextVertex(tracksExtra,massesExtra,vrtList2,*state);
2586 }
2587 }
2588 else { // m_JXSubVtx=false
2589 // Mother vertex includes just one subvertex (DisV) and JX tracks + extra track
2590 if(m_constrMainV) {
2591 vID3 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList2,*state,m_massMainV);
2592 } else {
2593 vID3 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList2,*state);
2594 }
2595 }
2596 if (m_constrJX && m_jxDaug_num>2) {
2597 std::vector<Trk::VertexID> cnstV;
2598 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
2599 ATH_MSG_WARNING("addMassConstraint for JX failed");
2600 }
2601 }
2602 if (m_constrJpsi) {
2603 std::vector<Trk::VertexID> cnstV;
2604 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
2605 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
2606 }
2607 }
2608 if (m_constrX && m_jxDaug_num==4) {
2609 std::vector<Trk::VertexID> cnstV;
2610 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksX,cnstV,*state,m_massX).isSuccess() ) {
2611 ATH_MSG_WARNING("addMassConstraint for X failed");
2612 }
2613 }
2614 // Do the work
2615 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) );
2616
2617 if (fit_result) {
2618 for(auto& v : fit_result->vertices()) {
2619 if(v->nTrackParticles()==0) {
2620 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
2621 v->setTrackParticleLinks(nullLinkVector);
2622 }
2623 }
2624 // reset links to original tracks
2625 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
2626
2627 // necessary to prevent memory leak
2628 fit_result->setSVOwnership(true);
2629
2630 // Chi2/DOF cut
2631 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
2632 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
2633
2634 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
2635 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
2636 size_t iMoth = cascadeVertices.size()-1;
2637 double lxy_SV1_sub = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
2638 double lxy_SV1 = m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[iMoth]);
2639
2640 if(chi2CutPassed && lxy_SV1>m_lxyDisV_cut && lxy_SV1_sub>m_lxyV0_cut) {
2641 chi2_V0_decor(*cascadeVertices[0]) = disVtx.V0vtx->chiSquared();
2642 ndof_V0_decor(*cascadeVertices[0]) = disVtx.V0vtx->numberDoF();
2643 if(disVtx.V0type==LAMBDA) {
2644 type_V0_decor(*cascadeVertices[0]) = "Lambda";
2645 }
2646 else if(disVtx.V0type==LAMBDABAR) {
2647 type_V0_decor(*cascadeVertices[0]) = "Lambdabar";
2648 }
2649 else if(disVtx.V0type==KS) {
2650 type_V0_decor(*cascadeVertices[0]) = "Ks";
2651 }
2652 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*disVtx.V0vtx) ? mAcc_gfit(*disVtx.V0vtx) : 0;
2653 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*disVtx.V0vtx) ? mAcc_gmass(*disVtx.V0vtx) : -1;
2654 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*disVtx.V0vtx) ? mAcc_gmasserr(*disVtx.V0vtx) : -1;
2655 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*disVtx.V0vtx) ? mAcc_gchisq(*disVtx.V0vtx) : 999999;
2656 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*disVtx.V0vtx) ? mAcc_gndof(*disVtx.V0vtx) : 0;
2657 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*disVtx.V0vtx) ? mAcc_gprob(*disVtx.V0vtx) : -1;
2658 trk_px.clear(); trk_py.clear(); trk_pz.clear();
2659 trk_px.push_back( disVtx.p4_V0track1.Px() ); trk_px.push_back( disVtx.p4_V0track2.Px() );
2660 trk_py.push_back( disVtx.p4_V0track1.Py() ); trk_py.push_back( disVtx.p4_V0track2.Py() );
2661 trk_pz.push_back( disVtx.p4_V0track1.Pz() ); trk_pz.push_back( disVtx.p4_V0track2.Pz() );
2662 trk_pxDeco(*cascadeVertices[0]) = trk_px;
2663 trk_pyDeco(*cascadeVertices[0]) = trk_py;
2664 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
2665 trk_px_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Px();
2666 trk_py_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Py();
2667 trk_pz_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Pz();
2668
2669 result.push_back( std::make_pair(fit_result.release(),nullptr) );
2670 }
2671 }
2672 } // loop over trackContainer
2673 } // m_extraTrk1MassHypo>0
2674
2675 if(pv_xAOD) {
2676 for(auto cascade_info_pair : result) {
2677 if(cascade_info_pair.first && cascade_info_pair.first->getParticleMoms().size()>0) {
2678 size_t index = cascade_info_pair.first->getParticleMoms().size() - 1;
2679 const std::vector<TLorentzVector> &mom = cascade_info_pair.first->getParticleMoms()[index];
2680 const Amg::MatrixX &cov = cascade_info_pair.first->getCovariance()[index];
2681 const xAOD::Vertex* mainVertex = cascade_info_pair.first->vertices()[index];
2682 xAOD::BPhysHypoHelper vtx(m_hypoName, mainVertex);
2683 bool isInDefaultPVCont = false;
2684 for(const xAOD::Vertex* pvVtx : *defaultPVContainer) {
2685 if(pv_xAOD == pvVtx) { isInDefaultPVCont = true; break; }
2686 }
2687 if(isInDefaultPVCont) vtx.setPv( pv_xAOD, defaultPVContainer, pvtype );
2688 else vtx.setPv( pv_xAOD, pvContainer, pvtype );
2689 if(origPv_xAOD) vtx.setOrigPv( origPv_xAOD, defaultPVContainer, pvtype );
2690 vtx.setLxy ( m_CascadeTools->lxy (mom, vtx.vtx(), pv_xAOD), pvtype );
2691 vtx.setLxyErr ( m_CascadeTools->lxyError (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2692 vtx.setA0 ( m_CascadeTools->a0 (mom, vtx.vtx(), pv_xAOD), pvtype );
2693 vtx.setA0Err ( m_CascadeTools->a0Error (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2694 vtx.setA0xy ( m_CascadeTools->a0xy (mom, vtx.vtx(), pv_xAOD), pvtype );
2695 vtx.setA0xyErr( m_CascadeTools->a0xyError(mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2696 vtx.setZ0 ( m_CascadeTools->a0z (mom, vtx.vtx(), pv_xAOD), pvtype );
2697 vtx.setZ0Err ( m_CascadeTools->a0zError (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2698 vtx.setRefitPVStatus( 0, pvtype );
2699 // Proper decay times
2700 vtx.setTau( m_CascadeTools->tau(mom, vtx.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
2701 vtx.setTauErr( m_CascadeTools->tauError(mom, cov, vtx.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
2702 vtx.setTau( m_CascadeTools->tau(mom, vtx.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
2703 vtx.setTauErr( m_CascadeTools->tauError(mom, cov, vtx.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
2704 }
2705 }
2706 }
2707
2708 return result;
2709 }
2710
2711 void JpsiXPlusDisplaced::fitV0Container(xAOD::VertexContainer* V0ContainerNew, const std::vector<const xAOD::TrackParticle*>& selectedTracks, const std::vector<const xAOD::TrackParticleContainer*>& trackCols) const {
2712 const EventContext& ctx = Gaudi::Hive::currentContext();
2713
2714 SG::AuxElement::Decorator<std::string> mDec_type("Type_V0Vtx");
2715 SG::AuxElement::Decorator<int> mDec_gfit("gamma_fit");
2716 SG::AuxElement::Decorator<float> mDec_gmass("gamma_mass");
2717 SG::AuxElement::Decorator<float> mDec_gmasserr("gamma_massError");
2718 SG::AuxElement::Decorator<float> mDec_gchisq("gamma_chisq");
2719 SG::AuxElement::Decorator<int> mDec_gndof("gamma_ndof");
2720 SG::AuxElement::Decorator<float> mDec_gprob("gamma_probability");
2721
2722 std::vector<const xAOD::TrackParticle*> posTracks;
2723 std::vector<const xAOD::TrackParticle*> negTracks;
2724 for(const xAOD::TrackParticle* TP : selectedTracks) {
2725 if(TP->charge()>0) posTracks.push_back(TP);
2726 else negTracks.push_back(TP);
2727 }
2728
2729 for(const xAOD::TrackParticle* TP1 : posTracks) {
2730 const Trk::Perigee& aPerigee1 = TP1->perigeeParameters();
2731 for(const xAOD::TrackParticle* TP2 : negTracks) {
2732 const Trk::Perigee& aPerigee2 = TP2->perigeeParameters();
2733 int sflag(0), errorcode(0);
2734 Amg::Vector3D startingPoint = m_vertexEstimator->getCirclesIntersectionPoint(&aPerigee1,&aPerigee2,sflag,errorcode);
2735 if (errorcode != 0) {startingPoint(0) = 0.0; startingPoint(1) = 0.0; startingPoint(2) = 0.0;}
2736
2737 if (errorcode == 0 || errorcode == 5 || errorcode == 6 || errorcode == 8) {
2738 Trk::PerigeeSurface perigeeSurface(startingPoint);
2739 const Trk::TrackParameters* extrapolatedPerigee1 = m_extrapolator->extrapolate(ctx,TP1->perigeeParameters(), perigeeSurface).release();
2740 const Trk::TrackParameters* extrapolatedPerigee2 = m_extrapolator->extrapolate(ctx,TP2->perigeeParameters(), perigeeSurface).release();
2741 std::vector<std::unique_ptr<const Trk::TrackParameters> > cleanup;
2742 if(!extrapolatedPerigee1) extrapolatedPerigee1 = &TP1->perigeeParameters();
2743 else cleanup.push_back(std::unique_ptr<const Trk::TrackParameters>(extrapolatedPerigee1));
2744 if(!extrapolatedPerigee2) extrapolatedPerigee2 = &TP2->perigeeParameters();
2745 else cleanup.push_back(std::unique_ptr<const Trk::TrackParameters>(extrapolatedPerigee2));
2746 if(extrapolatedPerigee1 && extrapolatedPerigee2) {
2747 bool pass = false;
2748 TLorentzVector v1; TLorentzVector v2;
2749 if(!pass) {
2750 v1.SetXYZM(extrapolatedPerigee1->momentum().x(),extrapolatedPerigee1->momentum().y(),extrapolatedPerigee1->momentum().z(),m_mass_proton);
2751 v2.SetXYZM(extrapolatedPerigee2->momentum().x(),extrapolatedPerigee2->momentum().y(),extrapolatedPerigee2->momentum().z(),m_mass_pion);
2752 if((v1+v2).M()>1030.0 && (v1+v2).M()<1200.0) pass = true;
2753 }
2754 if(!pass) {
2755 v1.SetXYZM(extrapolatedPerigee1->momentum().x(),extrapolatedPerigee1->momentum().y(),extrapolatedPerigee1->momentum().z(),m_mass_pion);
2756 v2.SetXYZM(extrapolatedPerigee2->momentum().x(),extrapolatedPerigee2->momentum().y(),extrapolatedPerigee2->momentum().z(),m_mass_proton);
2757 if((v1+v2).M()>1030.0 && (v1+v2).M()<1200.0) pass = true;
2758 }
2759 if(!pass) {
2760 v1.SetXYZM(extrapolatedPerigee1->momentum().x(),extrapolatedPerigee1->momentum().y(),extrapolatedPerigee1->momentum().z(),m_mass_pion);
2761 v2.SetXYZM(extrapolatedPerigee2->momentum().x(),extrapolatedPerigee2->momentum().y(),extrapolatedPerigee2->momentum().z(),m_mass_pion);
2762 if((v1+v2).M()>430.0 && (v1+v2).M()<565.0) pass = true;
2763 }
2764 if(pass) {
2765 std::vector<const xAOD::TrackParticle*> tracksV0;
2766 tracksV0.push_back(TP1); tracksV0.push_back(TP2);
2767 std::unique_ptr<xAOD::Vertex> V0vtx = std::unique_ptr<xAOD::Vertex>( m_iV0Fitter->fit(tracksV0, startingPoint) );
2768 if(V0vtx && V0vtx->chiSquared()>=0) {
2769 double chi2DOF = V0vtx->chiSquared()/V0vtx->numberDoF();
2770 if(chi2DOF>m_chi2cut_V0) continue;
2771
2772 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);
2773 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);
2774 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);
2775 if(massSig_V0_Lambda1<=massSig_V0_Lambda2 && massSig_V0_Lambda1<=massSig_V0_Ks) {
2776 mDec_type(*V0vtx.get()) = "Lambda";
2777 }
2778 else if(massSig_V0_Lambda2<=massSig_V0_Lambda1 && massSig_V0_Lambda2<=massSig_V0_Ks) {
2779 mDec_type(*V0vtx.get()) = "Lambdabar";
2780 }
2781 else if(massSig_V0_Ks<=massSig_V0_Lambda1 && massSig_V0_Ks<=massSig_V0_Lambda2) {
2782 mDec_type(*V0vtx.get()) = "Ks";
2783 }
2784
2785 int gamma_fit = 0; int gamma_ndof = 0; double gamma_chisq = 999999.;
2786 double gamma_prob = -1., gamma_mass = -1., gamma_massErr = -1.;
2787 std::unique_ptr<xAOD::Vertex> gammaVtx = std::unique_ptr<xAOD::Vertex>( m_iGammaFitter->fit(tracksV0, m_V0Tools->vtx(V0vtx.get())) );
2788 if (gammaVtx) {
2789 gamma_fit = 1;
2790 gamma_mass = m_V0Tools->invariantMass(gammaVtx.get(),m_mass_e,m_mass_e);
2791 gamma_massErr = m_V0Tools->invariantMassError(gammaVtx.get(),m_mass_e,m_mass_e);
2792 gamma_chisq = m_V0Tools->chisq(gammaVtx.get());
2793 gamma_ndof = m_V0Tools->ndof(gammaVtx.get());
2794 gamma_prob = m_V0Tools->vertexProbability(gammaVtx.get());
2795 }
2796 mDec_gfit(*V0vtx.get()) = gamma_fit;
2797 mDec_gmass(*V0vtx.get()) = gamma_mass;
2798 mDec_gmasserr(*V0vtx.get()) = gamma_massErr;
2799 mDec_gchisq(*V0vtx.get()) = gamma_chisq;
2800 mDec_gndof(*V0vtx.get()) = gamma_ndof;
2801 mDec_gprob(*V0vtx.get()) = gamma_prob;
2802
2803 xAOD::BPhysHelper V0_helper(V0vtx.get());
2804 V0_helper.setRefTrks(); // AOD only method
2805
2806 if(not trackCols.empty()){
2807 try {
2808 JpsiUpsilonCommon::RelinkVertexTracks(trackCols, V0vtx.get());
2809 } catch (std::runtime_error const& e) {
2810 ATH_MSG_ERROR(e.what());
2811 return;
2812 }
2813 }
2814
2815 V0ContainerNew->push_back(std::move(V0vtx));
2816 }
2817 }
2818 }
2819 }
2820 }
2821 }
2822 }
2823
2824 template<size_t NTracks>
2825 const xAOD::Vertex* JpsiXPlusDisplaced::FindVertex(const xAOD::VertexContainer* cont, const xAOD::Vertex* v) const {
2826 for (const xAOD::Vertex* v1 : *cont) {
2827 assert(v1->nTrackParticles() == NTracks);
2828 std::array<const xAOD::TrackParticle*, NTracks> a1;
2829 std::array<const xAOD::TrackParticle*, NTracks> a2;
2830 for(size_t i=0; i<NTracks; i++){
2831 a1[i] = v1->trackParticle(i);
2832 a2[i] = v->trackParticle(i);
2833 }
2834 std::sort(a1.begin(), a1.end());
2835 std::sort(a2.begin(), a2.end());
2836 if(a1 == a2) return v1;
2837 }
2838 return nullptr;
2839 }
2840}
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
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(const std::vector< const xAOD::TrackParticleContainer * > &trkcols, xAOD::Vertex *vtx)
Definition JpsiUpsilonCommon.cxx:126
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::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::getXiCandidate
XiCandidate getXiCandidate(const xAOD::Vertex *V0vtx, const V0Enum V0, const xAOD::TrackParticle *track3) const
Definition JpsiXPlusDisplaced.cxx:1179
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 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:1315
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::d0Pass
bool d0Pass(const xAOD::TrackParticle *track, const xAOD::Vertex *PV) const
Definition JpsiXPlusDisplaced.cxx:1168
DerivationFramework::JpsiXPlusDisplaced::FindVertex
const xAOD::Vertex * FindVertex(const xAOD::VertexContainer *cont, const xAOD::Vertex *v) const
Definition JpsiXPlusDisplaced.cxx:2825
DerivationFramework::JpsiXPlusDisplaced::m_v0VtxOutputKey
SG::WriteHandleKey< xAOD::VertexContainer > m_v0VtxOutputKey
Definition JpsiXPlusDisplaced.h:90
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_CascadeTools
ToolHandle< DerivationFramework::CascadeTools > m_CascadeTools
Definition JpsiXPlusDisplaced.h:192
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::getD0Candidate
MesonCandidate getD0Candidate(const xAOD::Vertex *JXvtx, const xAOD::TrackParticle *extraTrk1, const xAOD::TrackParticle *extraTrk2) const
Definition JpsiXPlusDisplaced.cxx:1284
DerivationFramework::JpsiXPlusDisplaced::fitV0Container
void fitV0Container(xAOD::VertexContainer *V0ContainerNew, const std::vector< const xAOD::TrackParticle * > &selectedTracks, const std::vector< const xAOD::TrackParticleContainer * > &trackCols) const
Definition JpsiXPlusDisplaced.cxx:2711
DerivationFramework::JpsiXPlusDisplaced::fitTracks
std::unique_ptr< xAOD::Vertex > fitTracks(const xAOD::TrackParticle *track1, const xAOD::TrackParticle *track2, const xAOD::TrackParticle *track3=nullptr) const
Definition JpsiXPlusDisplaced.cxx:1229
DerivationFramework::JpsiXPlusDisplaced::getDpmCandidate
MesonCandidate getDpmCandidate(const xAOD::Vertex *JXvtx, const xAOD::TrackParticle *extraTrk1, const xAOD::TrackParticle *extraTrk2, const xAOD::TrackParticle *extraTrk3) const
Definition JpsiXPlusDisplaced.cxx:1248
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_pvRefitter
ToolHandle< Analysis::PrimaryVertexRefitter > m_pvRefitter
Definition JpsiXPlusDisplaced.h:187
DerivationFramework::JpsiXPlusDisplaced::m_vertexEstimator
ToolHandle< InDet::VertexPointEstimator > m_vertexEstimator
Definition JpsiXPlusDisplaced.h:193
SG::AuxElement::Decorator
SG::Decorator< T, ALLOC > Decorator
Definition AuxElement.h:575
SG::AuxElement::Accessor
SG::Accessor< T, ALLOC > Accessor
Definition AuxElement.h:572
SG::ConstAccessor::isAvailable
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
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:415
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:270
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:249
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.
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
result
Generated on for ATLAS Offline Software by doxygen 1.14.0