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