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::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::Accessor<std::string> mAcc_type("Type_V0Vtx");
654 SG::Accessor<int> mAcc_gfit("gamma_fit");
655 SG::Accessor<float> mAcc_gmass("gamma_mass");
656 SG::Accessor<float> mAcc_gchisq("gamma_chisq");
657 SG::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(),
741 [](std::pair<const xAOD::Vertex*,V0Enum>& a, std::pair<const xAOD::Vertex*,V0Enum>& b)
742 { return a.first->chiSquared()/a.first->numberDoF() < b.first->chiSquared()/b.first->numberDoF(); } );
743 //
744 if(m_maxV0Candidates>0 && selectedV0Candidates.size()>m_maxV0Candidates) {
745 selectedV0Candidates.erase(selectedV0Candidates.begin()+m_maxV0Candidates, selectedV0Candidates.end());
746 }
747 if(selectedV0Candidates.size()==0) return StatusCode::SUCCESS;
748
749 std::vector<std::pair<Trk::VxCascadeInfo*, Trk::VxCascadeInfo*> > cascadeinfoContainer;
750 ATH_CHECK( performSearch(cascadeinfoContainer, selectedV0Candidates, tracksDisplaced, ctx) );
751
752 // sort and chop the main candidates
753 std::sort( cascadeinfoContainer.begin(), cascadeinfoContainer.end(),
754 [](std::pair<Trk::VxCascadeInfo*, Trk::VxCascadeInfo*> a, std::pair<Trk::VxCascadeInfo*, Trk::VxCascadeInfo*> b)
755 { return a.first->fitChi2()/a.first->nDoF() < b.first->fitChi2()/b.first->nDoF(); } );
756 if(m_maxMainVCandidates>0 && cascadeinfoContainer.size()>m_maxMainVCandidates) {
757 for(auto it=cascadeinfoContainer.begin()+m_maxMainVCandidates; it!=cascadeinfoContainer.end(); it++) {
758 delete it->first;
759 delete it->second;
760 }
761 cascadeinfoContainer.erase(cascadeinfoContainer.begin()+m_maxMainVCandidates, cascadeinfoContainer.end());
762 }
763 if(cascadeinfoContainer.size()==0) return StatusCode::SUCCESS;
764
765 SG::ReadHandle<xAOD::EventInfo> evt(m_eventInfo_key, ctx);
766 ATH_CHECK( evt.isValid() );
767 BPhysPVCascadeTools helper(&(*m_CascadeTools), evt.cptr());
768 helper.SetMinNTracksInPV(m_PV_minNTracks);
769
770 // Decorators for the cascade vertices
771 SG::Decorator<VertexLinkVector> CascadeLinksDecor("CascadeVertexLinks");
772 SG::Decorator<VertexLinkVector> PrecedingLinksDecor("PrecedingVertexLinks");
773 SG::Decorator<float> chi2_decor("ChiSquared");
774 SG::Decorator<int> ndof_decor("nDoF");
775 SG::Decorator<float> Pt_decor("Pt");
776 SG::Decorator<float> PtErr_decor("PtErr");
777
778 SG::Decorator<float> lxy_SV0_decor("lxy_SV0");
779 SG::Decorator<float> lxyErr_SV0_decor("lxyErr_SV0");
780 SG::Decorator<float> a0xy_SV0_decor("a0xy_SV0");
781 SG::Decorator<float> a0xyErr_SV0_decor("a0xyErr_SV0");
782 SG::Decorator<float> a0z_SV0_decor("a0z_SV0");
783 SG::Decorator<float> a0zErr_SV0_decor("a0zErr_SV0");
784
785 SG::Decorator<float> lxy_SV1_decor("lxy_SV1");
786 SG::Decorator<float> lxyErr_SV1_decor("lxyErr_SV1");
787 SG::Decorator<float> a0xy_SV1_decor("a0xy_SV1");
788 SG::Decorator<float> a0xyErr_SV1_decor("a0xyErr_SV1");
789 SG::Decorator<float> a0z_SV1_decor("a0z_SV1");
790 SG::Decorator<float> a0zErr_SV1_decor("a0zErr_SV1");
791
792 SG::Decorator<float> lxy_SV2_decor("lxy_SV2");
793 SG::Decorator<float> lxyErr_SV2_decor("lxyErr_SV2");
794 SG::Decorator<float> a0xy_SV2_decor("a0xy_SV2");
795 SG::Decorator<float> a0xyErr_SV2_decor("a0xyErr_SV2");
796 SG::Decorator<float> a0z_SV2_decor("a0z_SV2");
797 SG::Decorator<float> a0zErr_SV2_decor("a0zErr_SV2");
798
799 SG::Decorator<float> chi2_V2_decor("ChiSquared_V2");
800 SG::Decorator<int> ndof_V2_decor("nDoF_V2");
801
802 for(auto cascade_info_pair : cascadeinfoContainer) {
803 if(cascade_info_pair.first==nullptr) {
804 ATH_MSG_ERROR("CascadeInfo is null");
805 continue;
806 }
807
808 const std::vector<xAOD::Vertex*> &cascadeVertices = cascade_info_pair.first->vertices();
809 if(cascadeVertices.size() != topoN) ATH_MSG_ERROR("Incorrect number of vertices");
810 for(size_t i=0; i<topoN; i++) {
811 if(cascadeVertices[i]==nullptr) ATH_MSG_ERROR("Error null vertex");
812 }
813
814 cascade_info_pair.first->setSVOwnership(false); // Prevent Container from deleting vertices
815 const auto mainVertex = cascadeVertices[topoN-1]; // this is the mother vertex
816 const std::vector< std::vector<TLorentzVector> > &moms = cascade_info_pair.first->getParticleMoms();
817
818 // Identify the input JX
819 int ijx = m_JXSubVtx ? topoN-2 : topoN-1;
820 if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0) {
821 if(m_JXV0SubVtx) ijx = 1;
822 else ijx = topoN-1;
823 }
824 else if(m_extraTrk1MassHypo>0 && m_disVDaug_num==2) {
825 if(m_JXV0SubVtx || m_JXSubVtx) ijx = 1;
826 else ijx = topoN-1;
827 }
828 const xAOD::Vertex* jxVtx(nullptr);
829 if(m_jxDaug_num==4) jxVtx = FindVertex<4>(jxContainer.ptr(), cascadeVertices[ijx]);
830 else if(m_jxDaug_num==3) jxVtx = FindVertex<3>(jxContainer.ptr(), cascadeVertices[ijx]);
831 else jxVtx = FindVertex<2>(jxContainer.ptr(), cascadeVertices[ijx]);
832
833 xAOD::BPhysHypoHelper vtx(m_hypoName, mainVertex);
834
835 // Get refitted track momenta from all vertices, charged tracks only
836 BPhysPVCascadeTools::SetVectorInfo(vtx, cascade_info_pair.first);
837 vtx.setPass(true);
838
839 //
840 // Decorate main vertex
841 //
842 // mass, mass error
843 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/VxCascadeInfo.h
844 BPHYS_CHECK( vtx.setMass(m_CascadeTools->invariantMass(moms[topoN-1])) );
845 BPHYS_CHECK( vtx.setMassErr(m_CascadeTools->invariantMassError(moms[topoN-1],cascade_info_pair.first->getCovariance()[topoN-1])) );
846 // pt and pT error (the default pt of mainVertex is != the pt of the full cascade fit!)
847 Pt_decor(*mainVertex) = m_CascadeTools->pT(moms[topoN-1]);
848 PtErr_decor(*mainVertex) = m_CascadeTools->pTError(moms[topoN-1],cascade_info_pair.first->getCovariance()[topoN-1]);
849 // chi2 and ndof (the default chi2 of mainVertex is != the chi2 of the full cascade fit!)
850 chi2_decor(*mainVertex) = cascade_info_pair.first->fitChi2();
851 ndof_decor(*mainVertex) = cascade_info_pair.first->nDoF();
852
853 if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0) { // special cases
854 if(m_JXV0SubVtx) {
855 lxy_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
856 lxyErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->lxyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
857 a0z_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0z(moms[0],cascadeVertices[0],cascadeVertices[1]);
858 a0zErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0zError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
859 a0xy_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0xy(moms[0],cascadeVertices[0],cascadeVertices[1]);
860 a0xyErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0xyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
861 lxy_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->lxy(moms[1],cascadeVertices[1],mainVertex);
862 lxyErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->lxyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
863 a0z_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0z(moms[1],cascadeVertices[1],mainVertex);
864 a0zErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0zError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
865 a0xy_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0xy(moms[1],cascadeVertices[1],mainVertex);
866 a0xyErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0xyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
867 lxy_SV2_decor(*cascadeVertices[2]) = m_CascadeTools->lxy(moms[2],cascadeVertices[2],mainVertex);
868 lxyErr_SV2_decor(*cascadeVertices[2]) = m_CascadeTools->lxyError(moms[2],cascade_info_pair.first->getCovariance()[2],cascadeVertices[2],mainVertex);
869 a0z_SV2_decor(*cascadeVertices[2]) = m_CascadeTools->a0z(moms[2],cascadeVertices[2],mainVertex);
870 a0zErr_SV2_decor(*cascadeVertices[2]) = m_CascadeTools->a0zError(moms[2],cascade_info_pair.first->getCovariance()[2],cascadeVertices[2],mainVertex);
871 a0xy_SV2_decor(*cascadeVertices[2]) = m_CascadeTools->a0xy(moms[2],cascadeVertices[2],mainVertex);
872 a0xyErr_SV2_decor(*cascadeVertices[2]) = m_CascadeTools->a0xyError(moms[2],cascade_info_pair.first->getCovariance()[2],cascadeVertices[2],mainVertex);
873 }
874 else {
875 lxy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],mainVertex);
876 lxyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],mainVertex);
877 a0z_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0z(moms[0],cascadeVertices[0],mainVertex);
878 a0zErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0zError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],mainVertex);
879 a0xy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xy(moms[0],cascadeVertices[0],mainVertex);
880 a0xyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],mainVertex);
881 lxy_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->lxy(moms[1],cascadeVertices[1],mainVertex);
882 lxyErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->lxyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
883 a0z_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0z(moms[1],cascadeVertices[1],mainVertex);
884 a0zErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0zError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
885 a0xy_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0xy(moms[1],cascadeVertices[1],mainVertex);
886 a0xyErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0xyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
887 }
888 }
889 else if(m_extraTrk1MassHypo>0 && m_disVDaug_num==2) { // special cases
890 if(m_JXV0SubVtx) {
891 lxy_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
892 lxyErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->lxyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
893 a0z_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0z(moms[0],cascadeVertices[0],cascadeVertices[1]);
894 a0zErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0zError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
895 a0xy_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0xy(moms[0],cascadeVertices[0],cascadeVertices[1]);
896 a0xyErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0xyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
897 lxy_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->lxy(moms[1],cascadeVertices[1],mainVertex);
898 lxyErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->lxyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
899 a0z_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0z(moms[1],cascadeVertices[1],mainVertex);
900 a0zErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0zError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
901 a0xy_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0xy(moms[1],cascadeVertices[1],mainVertex);
902 a0xyErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0xyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
903 }
904 else {
905 lxy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],mainVertex);
906 lxyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],mainVertex);
907 a0z_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0z(moms[0],cascadeVertices[0],mainVertex);
908 a0zErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0zError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],mainVertex);
909 a0xy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xy(moms[0],cascadeVertices[0],mainVertex);
910 a0xyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],mainVertex);
911 if(m_JXSubVtx) {
912 lxy_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->lxy(moms[1],cascadeVertices[1],mainVertex);
913 lxyErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->lxyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
914 a0z_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0z(moms[1],cascadeVertices[1],mainVertex);
915 a0zErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0zError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
916 a0xy_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0xy(moms[1],cascadeVertices[1],mainVertex);
917 a0xyErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0xyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
918 }
919 }
920 }
921 else { // other normal cases
922 if(m_disVDaug_num==2) {
923 lxy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],mainVertex);
924 lxyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],mainVertex);
925 a0z_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0z(moms[0],cascadeVertices[0],mainVertex);
926 a0zErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0zError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],mainVertex);
927 a0xy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xy(moms[0],cascadeVertices[0],mainVertex);
928 a0xyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],mainVertex);
929 }
930 else {
931 lxy_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
932 lxyErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->lxyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
933 a0z_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0z(moms[0],cascadeVertices[0],cascadeVertices[1]);
934 a0zErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0zError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
935 a0xy_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0xy(moms[0],cascadeVertices[0],cascadeVertices[1]);
936 a0xyErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0xyError(moms[0],cascade_info_pair.first->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
937 lxy_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->lxy(moms[1],cascadeVertices[1],mainVertex);
938 lxyErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->lxyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
939 a0xy_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0z(moms[1],cascadeVertices[1],mainVertex);
940 a0xyErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0zError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
941 a0z_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0xy(moms[1],cascadeVertices[1],mainVertex);
942 a0zErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0xyError(moms[1],cascade_info_pair.first->getCovariance()[1],cascadeVertices[1],mainVertex);
943 }
944
945 if(m_JXSubVtx) {
946 lxy_SV2_decor(*cascadeVertices[ijx]) = m_CascadeTools->lxy(moms[ijx],cascadeVertices[ijx],mainVertex);
947 lxyErr_SV2_decor(*cascadeVertices[ijx]) = m_CascadeTools->lxyError(moms[ijx],cascade_info_pair.first->getCovariance()[ijx],cascadeVertices[ijx],mainVertex);
948 a0z_SV2_decor(*cascadeVertices[ijx]) = m_CascadeTools->a0z(moms[ijx],cascadeVertices[ijx],mainVertex);
949 a0zErr_SV2_decor(*cascadeVertices[ijx]) = m_CascadeTools->a0zError(moms[ijx],cascade_info_pair.first->getCovariance()[ijx],cascadeVertices[ijx],mainVertex);
950 a0xy_SV2_decor(*cascadeVertices[ijx]) = m_CascadeTools->a0xy(moms[ijx],cascadeVertices[ijx],mainVertex);
951 a0xyErr_SV2_decor(*cascadeVertices[ijx]) = m_CascadeTools->a0xyError(moms[ijx],cascade_info_pair.first->getCovariance()[ijx],cascadeVertices[ijx],mainVertex);
952 }
953 }
954
955 chi2_V2_decor(*cascadeVertices[ijx]) = m_V0Tools->chisq(jxVtx);
956 ndof_V2_decor(*cascadeVertices[ijx]) = m_V0Tools->ndof(jxVtx);
957
958 if(m_cascadeFitWithPV==0) {
959 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));
960 }
961
962 for(size_t i=0; i<topoN; i++) {
963 VtxWriteHandles[i].ptr()->push_back(cascadeVertices[i]);
964 }
965
966 // Set links to cascade vertices
967 VertexLinkVector cascadeVertexLinks;
968 VertexLink vertexLink1;
969 vertexLink1.setElement(cascadeVertices[0]);
970 vertexLink1.setStorableObject(*VtxWriteHandles[0].ptr());
971 if( vertexLink1.isValid() ) cascadeVertexLinks.push_back( vertexLink1 );
972 if(topoN>=3) {
973 VertexLink vertexLink2;
974 vertexLink2.setElement(cascadeVertices[1]);
975 vertexLink2.setStorableObject(*VtxWriteHandles[1].ptr());
976 if( vertexLink2.isValid() ) cascadeVertexLinks.push_back( vertexLink2 );
977 }
978 if(topoN==4) {
979 VertexLink vertexLink3;
980 vertexLink3.setElement(cascadeVertices[2]);
981 vertexLink3.setStorableObject(*VtxWriteHandles[2].ptr());
982 if( vertexLink3.isValid() ) cascadeVertexLinks.push_back( vertexLink3 );
983 }
984 CascadeLinksDecor(*mainVertex) = cascadeVertexLinks;
985
986 // process the posterior cascade with mass constraint for the main vertex
987 if(cascade_info_pair.second) {
988 const std::vector<xAOD::Vertex*> &cascadeVertices_mvc = cascade_info_pair.second->vertices();
989 if(cascadeVertices_mvc.size() != topoN) ATH_MSG_ERROR("Incorrect number of vertices (mvc)");
990 for(size_t i=0; i<topoN; i++) {
991 if(cascadeVertices_mvc[i]==nullptr) ATH_MSG_ERROR("Error null vertex (mvc)");
992 }
993 cascade_info_pair.second->setSVOwnership(false);
994 const auto mainVertex_mvc = cascadeVertices_mvc[topoN-1];
995 const std::vector< std::vector<TLorentzVector> > &moms_mvc = cascade_info_pair.second->getParticleMoms();
996 // Identify the input JX
997 int ijx_mvc = m_JXSubVtx ? topoN-2 : topoN-1;
998 if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0) {
999 if(m_JXV0SubVtx) ijx_mvc = 1;
1000 else ijx_mvc = topoN-1;
1001 }
1002 else if(m_extraTrk1MassHypo>0 && m_disVDaug_num==2) {
1003 if(m_JXV0SubVtx || m_JXSubVtx) ijx_mvc = 1;
1004 else ijx_mvc = topoN-1;
1005 }
1006 const xAOD::Vertex* jxVtx_mvc(nullptr);
1007 if(m_jxDaug_num==4) jxVtx_mvc = FindVertex<4>(jxContainer.ptr(), cascadeVertices_mvc[ijx_mvc]);
1008 else if(m_jxDaug_num==3) jxVtx_mvc = FindVertex<3>(jxContainer.ptr(), cascadeVertices_mvc[ijx_mvc]);
1009 else jxVtx_mvc = FindVertex<2>(jxContainer.ptr(), cascadeVertices_mvc[ijx_mvc]);
1010
1011 xAOD::BPhysHypoHelper vtx_mvc(m_hypoName, mainVertex_mvc);
1012 BPhysPVCascadeTools::SetVectorInfo(vtx_mvc, cascade_info_pair.second);
1013 vtx_mvc.setPass(true);
1014
1015 BPHYS_CHECK( vtx_mvc.setMass(m_CascadeTools->invariantMass(moms_mvc[topoN-1])) );
1016 BPHYS_CHECK( vtx_mvc.setMassErr(m_CascadeTools->invariantMassError(moms_mvc[topoN-1],cascade_info_pair.second->getCovariance()[topoN-1])) );
1017 Pt_decor(*mainVertex_mvc) = m_CascadeTools->pT(moms_mvc[topoN-1]);
1018 PtErr_decor(*mainVertex_mvc) = m_CascadeTools->pTError(moms_mvc[topoN-1],cascade_info_pair.second->getCovariance()[topoN-1]);
1019 chi2_decor(*mainVertex_mvc) = cascade_info_pair.second->fitChi2();
1020 ndof_decor(*mainVertex_mvc) = cascade_info_pair.second->nDoF();
1021
1022 if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0) { // special cases
1023 if(m_JXV0SubVtx) {
1024 lxy_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxy(moms_mvc[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1025 lxyErr_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1026 a0z_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0z(moms_mvc[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1027 a0zErr_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0zError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1028 a0xy_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xy(moms_mvc[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1029 a0xyErr_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1030 lxy_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1031 lxyErr_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1032 a0z_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0z(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1033 a0zErr_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0zError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1034 a0xy_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1035 a0xyErr_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1036 lxy_SV2_decor(*cascadeVertices_mvc[2]) = m_CascadeTools->lxy(moms_mvc[2],cascadeVertices_mvc[2],mainVertex_mvc);
1037 lxyErr_SV2_decor(*cascadeVertices_mvc[2]) = m_CascadeTools->lxyError(moms_mvc[2],cascade_info_pair.second->getCovariance()[2],cascadeVertices_mvc[2],mainVertex_mvc);
1038 a0z_SV2_decor(*cascadeVertices_mvc[2]) = m_CascadeTools->a0z(moms_mvc[2],cascadeVertices_mvc[2],mainVertex_mvc);
1039 a0zErr_SV2_decor(*cascadeVertices_mvc[2]) = m_CascadeTools->a0zError(moms_mvc[2],cascade_info_pair.second->getCovariance()[2],cascadeVertices_mvc[2],mainVertex_mvc);
1040 a0xy_SV2_decor(*cascadeVertices_mvc[2]) = m_CascadeTools->a0xy(moms_mvc[2],cascadeVertices_mvc[2],mainVertex_mvc);
1041 a0xyErr_SV2_decor(*cascadeVertices_mvc[2]) = m_CascadeTools->a0xyError(moms_mvc[2],cascade_info_pair.second->getCovariance()[2],cascadeVertices_mvc[2],mainVertex_mvc);
1042 }
1043 else {
1044 lxy_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxy(moms_mvc[0],cascadeVertices_mvc[0],mainVertex_mvc);
1045 lxyErr_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],mainVertex_mvc);
1046 a0z_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0z(moms_mvc[0],cascadeVertices_mvc[0],mainVertex_mvc);
1047 a0zErr_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0zError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],mainVertex_mvc);
1048 a0xy_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xy(moms_mvc[0],cascadeVertices_mvc[0],mainVertex_mvc);
1049 a0xyErr_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],mainVertex_mvc);
1050 lxy_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1051 lxyErr_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1052 a0z_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0z(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1053 a0zErr_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0zError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1054 a0xy_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1055 a0xyErr_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1056 }
1057 }
1058 else if(m_extraTrk1MassHypo>0 && m_disVDaug_num==2) { // special cases
1059 if(m_JXV0SubVtx) {
1060 lxy_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxy(moms_mvc[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1061 lxyErr_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1062 a0z_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0z(moms_mvc[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1063 a0zErr_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0zError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1064 a0xy_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xy(moms_mvc[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1065 a0xyErr_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1066 lxy_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1067 lxyErr_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1068 a0z_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0z(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1069 a0zErr_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0zError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1070 a0xy_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1071 a0xyErr_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1072 }
1073 else {
1074 lxy_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxy(moms_mvc[0],cascadeVertices_mvc[0],mainVertex_mvc);
1075 lxyErr_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],mainVertex_mvc);
1076 a0z_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0z(moms_mvc[0],cascadeVertices_mvc[0],mainVertex_mvc);
1077 a0zErr_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0zError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],mainVertex_mvc);
1078 a0xy_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xy(moms_mvc[0],cascadeVertices_mvc[0],mainVertex_mvc);
1079 a0xyErr_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],mainVertex_mvc);
1080 if(m_JXSubVtx) {
1081 lxy_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1082 lxyErr_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1083 a0z_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0z(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1084 a0zErr_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0zError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1085 a0xy_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1086 a0xyErr_SV2_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1087 }
1088 }
1089 }
1090 else { // other normal cases
1091 if(m_disVDaug_num==2) {
1092 lxy_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxy(moms_mvc[0],cascadeVertices_mvc[0],mainVertex_mvc);
1093 lxyErr_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],mainVertex_mvc);
1094 a0z_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0z(moms_mvc[0],cascadeVertices_mvc[0],mainVertex_mvc);
1095 a0zErr_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0zError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],mainVertex_mvc);
1096 a0xy_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xy(moms_mvc[0],cascadeVertices_mvc[0],mainVertex_mvc);
1097 a0xyErr_SV1_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],mainVertex_mvc);
1098 }
1099 else {
1100 lxy_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxy(moms_mvc[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1101 lxyErr_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->lxyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1102 a0z_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0z(moms_mvc[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1103 a0zErr_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0zError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1104 a0xy_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xy(moms_mvc[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1105 a0xyErr_SV0_decor(*cascadeVertices_mvc[0]) = m_CascadeTools->a0xyError(moms_mvc[0],cascade_info_pair.second->getCovariance()[0],cascadeVertices_mvc[0],cascadeVertices_mvc[1]);
1106 lxy_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1107 lxyErr_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->lxyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1108 a0xy_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0z(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1109 a0xyErr_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0zError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1110 a0z_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xy(moms_mvc[1],cascadeVertices_mvc[1],mainVertex_mvc);
1111 a0zErr_SV1_decor(*cascadeVertices_mvc[1]) = m_CascadeTools->a0xyError(moms_mvc[1],cascade_info_pair.second->getCovariance()[1],cascadeVertices_mvc[1],mainVertex_mvc);
1112 }
1113 if(m_JXSubVtx) {
1114 lxy_SV2_decor(*cascadeVertices_mvc[ijx_mvc]) = m_CascadeTools->lxy(moms_mvc[ijx_mvc],cascadeVertices_mvc[ijx_mvc],mainVertex_mvc);
1115 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);
1116 a0z_SV2_decor(*cascadeVertices_mvc[ijx_mvc]) = m_CascadeTools->a0z(moms_mvc[ijx_mvc],cascadeVertices_mvc[ijx_mvc],mainVertex_mvc);
1117 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);
1118 a0xy_SV2_decor(*cascadeVertices_mvc[ijx_mvc]) = m_CascadeTools->a0xy(moms_mvc[ijx_mvc],cascadeVertices_mvc[ijx_mvc],mainVertex_mvc);
1119 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);
1120 }
1121 }
1122 chi2_V2_decor(*cascadeVertices_mvc[ijx_mvc]) = m_V0Tools->chisq(jxVtx_mvc);
1123 ndof_V2_decor(*cascadeVertices_mvc[ijx_mvc]) = m_V0Tools->ndof(jxVtx_mvc);
1124
1125 if(m_cascadeFitWithPV==0) {
1126 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));
1127 }
1128
1129 for(size_t i=0; i<topoN; i++) {
1130 VtxWriteHandles_mvc[i].ptr()->push_back(cascadeVertices_mvc[i]);
1131 }
1132
1133 VertexLinkVector cascadeVertexLinks_mvc;
1134 VertexLink vertexLink1_mvc;
1135 vertexLink1_mvc.setElement(cascadeVertices_mvc[0]);
1136 vertexLink1_mvc.setStorableObject(*VtxWriteHandles_mvc[0].ptr());
1137 if( vertexLink1_mvc.isValid() ) cascadeVertexLinks_mvc.push_back( vertexLink1_mvc );
1138 if(topoN>=3) {
1139 VertexLink vertexLink2_mvc;
1140 vertexLink2_mvc.setElement(cascadeVertices_mvc[1]);
1141 vertexLink2_mvc.setStorableObject(*VtxWriteHandles_mvc[1].ptr());
1142 if( vertexLink2_mvc.isValid() ) cascadeVertexLinks_mvc.push_back( vertexLink2_mvc );
1143 }
1144 if(topoN==4) {
1145 VertexLink vertexLink3_mvc;
1146 vertexLink3_mvc.setElement(cascadeVertices_mvc[2]);
1147 vertexLink3_mvc.setStorableObject(*VtxWriteHandles_mvc[2].ptr());
1148 if( vertexLink3_mvc.isValid() ) cascadeVertexLinks_mvc.push_back( vertexLink3_mvc );
1149 }
1150 CascadeLinksDecor(*mainVertex_mvc) = cascadeVertexLinks_mvc;
1151
1152 VertexLinkVector precedingVertexLinks_mvc;
1153 VertexLink vertexLink_mvc;
1154 vertexLink_mvc.setElement(cascadeVertices[topoN-1]);
1155 vertexLink_mvc.setStorableObject(*VtxWriteHandles[topoN-1].ptr());
1156 if( vertexLink_mvc.isValid() ) precedingVertexLinks_mvc.push_back( vertexLink_mvc );
1157 PrecedingLinksDecor(*mainVertex_mvc) = precedingVertexLinks_mvc;
1158 } // cascade_info_pair.second!=0
1159 } // loop over cascadeinfoContainer
1160
1161 // Deleting cascadeinfo since this won't be stored.
1162 for (auto cascade_info_pair : cascadeinfoContainer) {
1163 delete cascade_info_pair.first;
1164 delete cascade_info_pair.second;
1165 }
1166
1167 return StatusCode::SUCCESS;
1168 }
1169
1170 bool JpsiXPlusDisplaced::d0Pass(const EventContext& ctx, const xAOD::TrackParticle* track, const xAOD::Vertex* PV) const {
1171 bool pass = false;
1172 std::unique_ptr<Trk::Perigee> per = m_trackToVertexTool->perigeeAtVertex(ctx, *track, PV->position());
1173 if(!per) return pass;
1174 double d0 = per->parameters()[Trk::d0];
1175 double sig_d0 = sqrt((*per->covariance())(0,0));
1176 if(std::abs(d0/sig_d0) > m_d0_cut) pass = true;
1177 return pass;
1178 }
1179
1180 JpsiXPlusDisplaced::XiCandidate JpsiXPlusDisplaced::getXiCandidate(const EventContext& ctx, const xAOD::Vertex* V0vtx, const V0Enum V0, const xAOD::TrackParticle* track3) const {
1181 XiCandidate disVtx;
1182
1183 std::vector<const xAOD::TrackParticle*> tracksV0;
1184 tracksV0.reserve(V0vtx->nTrackParticles());
1185 for(size_t i=0; i<V0vtx->nTrackParticles(); i++) tracksV0.push_back(V0vtx->trackParticle(i));
1186 std::vector<double> massesV0;
1187 if(V0==LAMBDA) massesV0 = m_massesV0_ppi;
1188 else if(V0==LAMBDABAR) massesV0 = m_massesV0_pip;
1189 else if(V0==KS) massesV0 = m_massesV0_pipi;
1190
1191 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
1192 int robustness = 0;
1193 m_iVertexFitter->setRobustness(robustness, *state);
1194 std::vector<Trk::VertexID> vrtList;
1195 // V0 vertex
1196 Trk::VertexID vID = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1197 vrtList.push_back(vID);
1198 // Mother vertex
1199 std::vector<const xAOD::TrackParticle*> tracksDis3{track3};
1200 std::vector<double> massesDis3{m_disVDaug3MassHypo};
1201 m_iVertexFitter->nextVertex(tracksDis3,massesDis3,vrtList,*state);
1202 // Do the work
1203 std::unique_ptr<Trk::VxCascadeInfo> cascade_info = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state) );
1204 if(cascade_info) {
1205 cascade_info->setSVOwnership(true);
1206 double chi2NDF = cascade_info->fitChi2()/cascade_info->nDoF();
1207 if(m_chi2cut_DisV<=0 || chi2NDF < m_chi2cut_DisV) {
1208 const std::vector<std::vector<TLorentzVector> > &moms = cascade_info->getParticleMoms();
1209 const std::vector<xAOD::Vertex*> &cascadeVertices = cascade_info->vertices();
1210 double lxy_SV1_sub = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
1211 if(lxy_SV1_sub > 0.2) {
1212 TLorentzVector totalMom;
1213 for(size_t it=0; it<moms[1].size(); it++) totalMom += moms[1][it];
1214 double disV_mass = totalMom.M();
1215 if(m_useImprovedMass) {
1216 if((V0==LAMBDA || V0==LAMBDABAR) && m_massLd>0) disV_mass += - moms[1][1].M() + m_massLd;
1217 else if(V0==KS && m_massKs>0) disV_mass += - moms[1][1].M() + m_massKs;
1218 }
1219 if(disV_mass>m_DisplacedMassLower && disV_mass<m_DisplacedMassUpper) {
1220 disVtx.track = track3; disVtx.V0vtx = V0vtx; disVtx.V0type = V0;
1221 disVtx.chi2NDF = chi2NDF; disVtx.p4_V0track1 = moms[0][0];
1222 disVtx.p4_V0track2 = moms[0][1]; disVtx.p4_disVtrack = moms[1][0];
1223 }
1224 }
1225 }
1226 }
1227 return disVtx;
1228 }
1229
1230 std::unique_ptr<xAOD::Vertex> JpsiXPlusDisplaced::fitTracks(const EventContext& ctx, const xAOD::TrackParticle* track1, const xAOD::TrackParticle* track2, const xAOD::TrackParticle* track3) const {
1231 // Starting point
1232 const Trk::Perigee& aPerigee1 = track1->perigeeParameters();
1233 const Trk::Perigee& aPerigee2 = track2->perigeeParameters();
1234 int sflag(0), errorcode(0);
1235 Amg::Vector3D startingPoint = m_vertexEstimator->getCirclesIntersectionPoint(&aPerigee1,&aPerigee2,sflag,errorcode);
1236 if(errorcode) startingPoint(0) = startingPoint(1) = startingPoint(2) = 0.0;
1237 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
1238 // do the fit
1239 if(track3) {
1240 return m_iVertexFitter->fit(std::vector<const xAOD::TrackParticle*>{track1,track2,track3}, startingPoint, *state);
1241 }
1242 else {
1243 return m_iVertexFitter->fit(std::vector<const xAOD::TrackParticle*>{track1,track2}, startingPoint, *state);
1244 }
1245 }
1246
1247 JpsiXPlusDisplaced::MesonCandidate JpsiXPlusDisplaced::getDpmCandidate(const EventContext& ctx, const xAOD::Vertex* JXvtx, const xAOD::TrackParticle* extraTrk1, const xAOD::TrackParticle* extraTrk2, const xAOD::TrackParticle* extraTrk3) const {
1248 MesonCandidate Dpm;
1249 // Check overlap
1250 std::vector<const xAOD::TrackParticle*> tracksJX;
1251 tracksJX.reserve(JXvtx->nTrackParticles());
1252 for(size_t i=0; i<JXvtx->nTrackParticles(); i++) tracksJX.push_back(JXvtx->trackParticle(i));
1253
1254 TLorentzVector tmp1, tmp2, tmp3;
1255 tmp1.SetPtEtaPhiM(extraTrk1->pt(),extraTrk1->eta(),extraTrk1->phi(),m_extraTrk1MassHypo);
1256 tmp2.SetPtEtaPhiM(extraTrk2->pt(),extraTrk2->eta(),extraTrk2->phi(),m_extraTrk2MassHypo);
1257 tmp3.SetPtEtaPhiM(extraTrk3->pt(),extraTrk3->eta(),extraTrk3->phi(),m_extraTrk3MassHypo);
1258 if((tmp1+tmp2+tmp3).M() < m_DpmMassLower || (tmp1+tmp2+tmp3).M() > m_DpmMassUpper) return Dpm;
1259
1260 std::unique_ptr<xAOD::Vertex> vtx = fitTracks(ctx, extraTrk1, extraTrk2, extraTrk3);
1261 if(vtx) {
1262 double chi2NDF = vtx->chiSquared()/vtx->numberDoF();
1263 if(m_chi2cut_Dpm<=0.0 || chi2NDF < m_chi2cut_Dpm) {
1264 double lxyDpm = m_V0Tools->lxy(vtx.get(),JXvtx);
1265 if(lxyDpm>m_lxyDpm_cut) {
1266 Dpm.extraTrack1 = extraTrk1; Dpm.extraTrack2 = extraTrk2; Dpm.extraTrack3 = extraTrk3;
1267 Dpm.chi2NDF = chi2NDF;
1268 TVector3 tot_pt; TVector3 tmp;
1269 for(size_t i=0; i<vtx->vxTrackAtVertex().size(); ++i) {
1270 const Trk::TrackParameters* aPerigee = vtx->vxTrackAtVertex()[i].perigeeAtVertex();
1271 if(aPerigee) {
1272 tmp.SetXYZ(aPerigee->momentum()[Trk::px],aPerigee->momentum()[Trk::py],aPerigee->momentum()[Trk::pz]);
1273 tot_pt += tmp;
1274 }
1275 }
1276 Dpm.pt = tot_pt.Pt();
1277 }
1278 }
1279 }
1280 return Dpm;
1281 }
1282
1283 JpsiXPlusDisplaced::MesonCandidate JpsiXPlusDisplaced::getD0Candidate(const EventContext& ctx, const xAOD::Vertex* JXvtx, const xAOD::TrackParticle* extraTrk1, const xAOD::TrackParticle* extraTrk2) const {
1284 MesonCandidate D0;
1285
1286 TLorentzVector tmp1, tmp2;
1287 tmp1.SetPtEtaPhiM(extraTrk1->pt(),extraTrk1->eta(),extraTrk1->phi(),m_extraTrk1MassHypo);
1288 tmp2.SetPtEtaPhiM(extraTrk2->pt(),extraTrk2->eta(),extraTrk2->phi(),m_extraTrk2MassHypo);
1289 if((tmp1+tmp2).M() < m_D0MassLower || (tmp1+tmp2).M() > m_D0MassUpper) return D0;
1290
1291 std::unique_ptr<xAOD::Vertex> vtx = fitTracks(ctx, extraTrk1, extraTrk2);
1292 if(vtx) {
1293 double chi2NDF = vtx->chiSquared()/vtx->numberDoF();
1294 if(m_chi2cut_D0<=0.0 || chi2NDF < m_chi2cut_D0) {
1295 double lxyD0 = m_V0Tools->lxy(vtx.get(),JXvtx);
1296 if(lxyD0>m_lxyD0_cut) {
1297 D0.extraTrack1 = extraTrk1; D0.extraTrack2 = extraTrk2;
1298 D0.chi2NDF = chi2NDF;
1299 TVector3 tot_pt; TVector3 tmp;
1300 for(size_t i=0; i<vtx->vxTrackAtVertex().size(); ++i) {
1301 const Trk::TrackParameters* aPerigee = vtx->vxTrackAtVertex()[i].perigeeAtVertex();
1302 if(aPerigee) {
1303 tmp.SetXYZ(aPerigee->momentum()[Trk::px],aPerigee->momentum()[Trk::py],aPerigee->momentum()[Trk::pz]);
1304 tot_pt += tmp;
1305 }
1306 }
1307 D0.pt = tot_pt.Pt();
1308 }
1309 }
1310 }
1311 return D0;
1312 }
1313
1314 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 {
1315 std::vector<std::pair<Trk::VxCascadeInfo*,Trk::VxCascadeInfo*> > result;
1316
1317 std::vector<const xAOD::TrackParticle*> tracksJX;
1318 tracksJX.reserve(JXvtx->nTrackParticles());
1319 for(size_t i=0; i<JXvtx->nTrackParticles(); i++) tracksJX.push_back(JXvtx->trackParticle(i));
1320 if (tracksJX.size() != massesJX.size()) {
1321 ATH_MSG_ERROR("Problems with JX input: number of tracks or track mass inputs is not correct!");
1322 return result;
1323 }
1324 // Check identical tracks in input
1325 if(std::find(tracksJX.cbegin(), tracksJX.cend(), V0vtx->trackParticle(0)) != tracksJX.cend()) return result;
1326 if(std::find(tracksJX.cbegin(), tracksJX.cend(), V0vtx->trackParticle(1)) != tracksJX.cend()) return result;
1327 std::vector<const xAOD::TrackParticle*> tracksV0;
1328 tracksV0.reserve(V0vtx->nTrackParticles());
1329 for(size_t j=0; j<V0vtx->nTrackParticles(); j++) tracksV0.push_back(V0vtx->trackParticle(j));
1330
1331 std::vector<const xAOD::TrackParticle*> tracksJpsi{tracksJX[0], tracksJX[1]};
1332 std::vector<const xAOD::TrackParticle*> tracksX;
1333 if(m_jxDaug_num>=3) tracksX.push_back(tracksJX[2]);
1334 if(m_jxDaug_num==4) tracksX.push_back(tracksJX[3]);
1335
1336 std::vector<double> massesV0;
1337 if(V0==LAMBDA) {
1338 massesV0 = m_massesV0_ppi;
1339 }
1340 else if(V0==LAMBDABAR) {
1341 massesV0 = m_massesV0_pip;
1342 }
1343 else if(V0==KS) {
1344 massesV0 = m_massesV0_pipi;
1345 }
1346
1347 TLorentzVector p4_moth, p4_v0, tmp;
1348 for(size_t it=0; it<JXvtx->nTrackParticles(); it++) {
1349 tmp.SetPtEtaPhiM(JXvtx->trackParticle(it)->pt(), JXvtx->trackParticle(it)->eta(), JXvtx->trackParticle(it)->phi(), massesJX[it]);
1350 p4_moth += tmp;
1351 }
1352 xAOD::BPhysHelper V0_helper(V0vtx);
1353 for(int it=0; it<V0_helper.nRefTrks(); it++) {
1354 p4_moth += V0_helper.refTrk(it,massesV0[it]);
1355 p4_v0 += V0_helper.refTrk(it,massesV0[it]);
1356 }
1357
1358 xAOD::BPhysHelper::pv_type pvtype = xAOD::BPhysHelper::PV_MIN_A0;
1359 if(m_cascadeFitWithPV==1) pvtype = xAOD::BPhysHelper::PV_MAX_SUM_PT2;
1360 else if(m_cascadeFitWithPV==2) pvtype = xAOD::BPhysHelper::PV_MIN_A0;
1361 else if(m_cascadeFitWithPV==3) pvtype = xAOD::BPhysHelper::PV_MIN_Z0;
1362 else if(m_cascadeFitWithPV==4) pvtype = xAOD::BPhysHelper::PV_MIN_Z0_BA;
1363 xAOD::BPhysHelper JX_helper(JXvtx);
1364 const xAOD::Vertex* origPv_xAOD = (m_cascadeFitWithPV>=1 && m_cascadeFitWithPV<=4) ? JX_helper.origPv(pvtype) : nullptr;
1365 const xAOD::Vertex* pv_xAOD = (m_cascadeFitWithPV>=1 && m_cascadeFitWithPV<=4) ? JX_helper.pv(pvtype) : nullptr;
1366 std::unique_ptr<Trk::RecVertex> pv_AOD;
1367 if(pv_xAOD) pv_AOD = std::make_unique<Trk::RecVertex>(pv_xAOD->position(),pv_xAOD->covariancePosition(),pv_xAOD->numberDoF(),pv_xAOD->chiSquared());
1368
1369 SG::Decorator<float> chi2_V1_decor("ChiSquared_V1");
1370 SG::Decorator<int> ndof_V1_decor("nDoF_V1");
1371 SG::Decorator<std::string> type_V1_decor("Type_V1");
1372
1373 SG::Accessor<int> mAcc_gfit("gamma_fit");
1374 SG::Accessor<float> mAcc_gmass("gamma_mass");
1375 SG::Accessor<float> mAcc_gmasserr("gamma_massError");
1376 SG::Accessor<float> mAcc_gchisq("gamma_chisq");
1377 SG::Accessor<int> mAcc_gndof("gamma_ndof");
1378 SG::Accessor<float> mAcc_gprob("gamma_probability");
1379
1380 SG::Decorator<int> mDec_gfit("gamma_fit");
1381 SG::Decorator<float> mDec_gmass("gamma_mass");
1382 SG::Decorator<float> mDec_gmasserr("gamma_massError");
1383 SG::Decorator<float> mDec_gchisq("gamma_chisq");
1384 SG::Decorator<int> mDec_gndof("gamma_ndof");
1385 SG::Decorator<float> mDec_gprob("gamma_probability");
1386 SG::Decorator< std::vector<float> > trk_pxDeco("TrackPx_V0nc");
1387 SG::Decorator< std::vector<float> > trk_pyDeco("TrackPy_V0nc");
1388 SG::Decorator< std::vector<float> > trk_pzDeco("TrackPz_V0nc");
1389
1390 std::vector<float> trk_px;
1391 std::vector<float> trk_py;
1392 std::vector<float> trk_pz;
1393
1394 if(m_extraTrk1MassHypo<=0) {
1395 double main_mass = p4_moth.M();
1396 if(m_useImprovedMass) {
1397 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_v0).M() + m_massJpsi;
1398 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_v0).M() + m_massJX;
1399 if((V0==LAMBDA || V0==LAMBDABAR) && m_massLd>0) main_mass += - p4_v0.M() + m_massLd;
1400 else if(V0==KS && m_massKs>0) main_mass += - p4_v0.M() + m_massKs;
1401 }
1402 if (main_mass < m_MassLower || main_mass > m_MassUpper) return result;
1403
1404 // Apply the user's settings to the fitter
1405 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
1406 // Robustness: http://cdsweb.cern.ch/record/685551
1407 int robustness = 0;
1408 m_iVertexFitter->setRobustness(robustness, *state);
1409 // Build up the topology
1410 // Vertex list
1411 std::vector<Trk::VertexID> vrtList;
1412 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
1413 // V0 vertex
1414 Trk::VertexID vID1;
1415 if (m_constrV0) {
1416 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,V0==KS?m_massKs:m_massLd);
1417 } else {
1418 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1419 }
1420 vrtList.push_back(vID1);
1421 Trk::VertexID vID2;
1422 if(m_JXSubVtx) {
1423 // JX vertex
1424 if (m_constrJX && m_jxDaug_num>2) {
1425 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state,m_massJX);
1426 } else {
1427 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state);
1428 }
1429 vrtList.push_back(vID2);
1430 // Mother vertex including JX and V0
1431 std::vector<const xAOD::TrackParticle*> tp;
1432 std::vector<double> tp_masses;
1433 if(m_constrMainV) {
1434 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList,*state,m_massMainV);
1435 } else {
1436 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList,*state);
1437 }
1438 }
1439 else { // m_JXSubVtx=false
1440 // Mother vertex including JX and V0
1441 if(m_constrMainV) {
1442 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state,m_massMainV);
1443 } else {
1444 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
1445 }
1446 if (m_constrJX && m_jxDaug_num>2) {
1447 std::vector<Trk::VertexID> cnstV;
1448 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
1449 ATH_MSG_WARNING("addMassConstraint for JX failed");
1450 }
1451 }
1452 }
1453 if (m_constrJpsi) {
1454 std::vector<Trk::VertexID> cnstV;
1455 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
1456 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
1457 }
1458 }
1459 if (m_constrX && m_jxDaug_num==4) {
1460 std::vector<Trk::VertexID> cnstV;
1461 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksX,cnstV,*state,m_massX).isSuccess() ) {
1462 ATH_MSG_WARNING("addMassConstraint for X failed");
1463 }
1464 }
1465 // Do the work
1466 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) );
1467
1468 if (fit_result) {
1469 for(auto& v : fit_result->vertices()) {
1470 if(v->nTrackParticles()==0) {
1471 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
1472 v->setTrackParticleLinks(nullLinkVector);
1473 }
1474 }
1475 // reset links to original tracks
1476 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
1477
1478 // necessary to prevent memory leak
1479 fit_result->setSVOwnership(true);
1480
1481 // Chi2/DOF cut
1482 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
1483 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
1484
1485 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
1486 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
1487 size_t iMoth = cascadeVertices.size()-1;
1488 double lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[iMoth]);
1489 if(chi2CutPassed && lxy_SV1>m_lxyV0_cut) {
1490 chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
1491 ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
1492 if(V0==LAMBDA) {
1493 type_V1_decor(*cascadeVertices[0]) = "Lambda";
1494 }
1495 else if(V0==LAMBDABAR) {
1496 type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
1497 }
1498 else if(V0==KS) {
1499 type_V1_decor(*cascadeVertices[0]) = "Ks";
1500 }
1501 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
1502 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
1503 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
1504 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
1505 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
1506 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
1507 trk_px.clear(); trk_py.clear(); trk_pz.clear();
1508 trk_px.reserve(V0_helper.nRefTrks());
1509 trk_py.reserve(V0_helper.nRefTrks());
1510 trk_pz.reserve(V0_helper.nRefTrks());
1511 for(auto&& vec3 : V0_helper.refTrks()) {
1512 trk_px.push_back( vec3.Px() );
1513 trk_py.push_back( vec3.Py() );
1514 trk_pz.push_back( vec3.Pz() );
1515 }
1516 trk_pxDeco(*cascadeVertices[0]) = std::move(trk_px);
1517 trk_pyDeco(*cascadeVertices[0]) = std::move(trk_py);
1518 trk_pzDeco(*cascadeVertices[0]) = std::move(trk_pz);
1519
1520 result.push_back( std::make_pair(fit_result.release(),nullptr) );
1521 }
1522 }
1523 } // for m_extraTrk1MassHypo<=0
1524 else if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo<=0) {
1525 std::vector<double> massesExtra{m_extraTrk1MassHypo};
1526 std::vector<double> massesJXExtra = massesJX; massesJXExtra.push_back(m_extraTrk1MassHypo);
1527
1528 for(const xAOD::TrackParticle* tpExtra : *trackContainer) {
1529 if( tpExtra->pt()<m_extraTrk1MinPt ) continue;
1530 if( !m_trkSelector->decision(*tpExtra, nullptr) ) continue;
1531 // Check identical tracks in input
1532 if(std::find(tracksJX.cbegin(),tracksJX.cend(),tpExtra) != tracksJX.cend()) continue;
1533 if(std::find(tracksV0.cbegin(),tracksV0.cend(),tpExtra) != tracksV0.cend()) continue;
1534
1535 TLorentzVector tmp;
1536 tmp.SetPtEtaPhiM(tpExtra->pt(),tpExtra->eta(),tpExtra->phi(),m_extraTrk1MassHypo);
1537 double main_mass = (p4_moth+tmp).M();
1538 if(m_useImprovedMass) {
1539 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_v0).M() + m_massJpsi;
1540 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_v0).M() + m_massJX;
1541 if((V0==LAMBDA || V0==LAMBDABAR) && m_massLd>0) main_mass += - p4_v0.M() + m_massLd;
1542 else if(V0==KS && m_massKs>0) main_mass += - p4_v0.M() + m_massKs;
1543 }
1544 if(main_mass < m_MassLower || main_mass > m_MassUpper) continue;
1545
1546 std::vector<const xAOD::TrackParticle*> tracksExtra{tpExtra};
1547 std::vector<const xAOD::TrackParticle*> tracksJXExtra = tracksJX; tracksJXExtra.push_back(tpExtra);
1548
1549 // Apply the user's settings to the fitter
1550 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
1551 // Robustness: http://cdsweb.cern.ch/record/685551
1552 int robustness = 0;
1553 m_iVertexFitter->setRobustness(robustness, *state);
1554 // Build up the topology
1555 // Vertex list
1556 std::vector<Trk::VertexID> vrtList;
1557 std::vector<Trk::VertexID> vrtList2;
1558 Trk::VertexID vID2;
1559 if(m_JXV0SubVtx) {
1560 // V0 vertex
1561 Trk::VertexID vID1;
1562 if (m_constrV0) {
1563 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,V0==KS?m_massKs:m_massLd);
1564 } else {
1565 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1566 }
1567 vrtList.push_back(vID1);
1568 // JX+V0 vertex
1569 if(m_constrJXV0) {
1570 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state,m_massJXV0);
1571 } else {
1572 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
1573 }
1574 vrtList2.push_back(vID2);
1575 // Main vertex
1576 if(m_constrMainV) {
1577 m_iVertexFitter->nextVertex(tracksExtra,massesExtra,vrtList2,*state,m_massMainV);
1578 } else {
1579 m_iVertexFitter->nextVertex(tracksExtra,massesExtra,vrtList2,*state);
1580 }
1581 }
1582 else { // m_JXV0SubVtx==false
1583 // V0 vertex
1584 Trk::VertexID vID1;
1585 if (m_constrV0) {
1586 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,V0==KS?m_massKs:m_massLd);
1587 } else {
1588 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1589 }
1590 vrtList.push_back(vID1);
1591 if(m_JXSubVtx) {
1592 // JX vertex
1593 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state);
1594 vrtList.push_back(vID2);
1595 // Mother vertex
1596 if(m_constrMainV) {
1597 m_iVertexFitter->nextVertex(tracksExtra,massesExtra,vrtList,*state,m_massMainV);
1598 } else {
1599 m_iVertexFitter->nextVertex(tracksExtra,massesExtra,vrtList,*state);
1600 }
1601 }
1602 else { // m_JXSubVtx=false
1603 // Mother vertex includes one subvertex (V0) and JX tracks + extra track
1604 if(m_constrMainV) {
1605 vID2 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList,*state,m_massMainV);
1606 } else {
1607 vID2 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList,*state);
1608 }
1609 }
1610 }
1611 if (m_constrJX && m_jxDaug_num>2) {
1612 std::vector<Trk::VertexID> cnstV;
1613 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
1614 ATH_MSG_WARNING("addMassConstraint for JX failed");
1615 }
1616 }
1617 if (m_constrJpsi) {
1618 std::vector<Trk::VertexID> cnstV;
1619 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
1620 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
1621 }
1622 }
1623 if (m_constrX && m_jxDaug_num==4) {
1624 std::vector<Trk::VertexID> cnstV;
1625 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksX,cnstV,*state,m_massX).isSuccess() ) {
1626 ATH_MSG_WARNING("addMassConstraint for X failed");
1627 }
1628 }
1629 // Do the work
1630 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) );
1631
1632 if (fit_result) {
1633 for(auto& v : fit_result->vertices()) {
1634 if(v->nTrackParticles()==0) {
1635 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
1636 v->setTrackParticleLinks(nullLinkVector);
1637 }
1638 }
1639 // reset links to original tracks
1640 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
1641
1642 // necessary to prevent memory leak
1643 fit_result->setSVOwnership(true);
1644
1645 // Chi2/DOF cut
1646 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
1647 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
1648 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
1649 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
1650 size_t iMoth = cascadeVertices.size()-1;
1651 double lxy_SV1(0);
1652 if(m_JXV0SubVtx) {
1653 lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
1654 }
1655 else {
1656 lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[iMoth]);
1657 }
1658 if(chi2CutPassed && lxy_SV1>m_lxyV0_cut) {
1659 chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
1660 ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
1661 if(V0==LAMBDA) {
1662 type_V1_decor(*cascadeVertices[0]) = "Lambda";
1663 }
1664 else if(V0==LAMBDABAR) {
1665 type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
1666 }
1667 else if(V0==KS) {
1668 type_V1_decor(*cascadeVertices[0]) = "Ks";
1669 }
1670 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
1671 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
1672 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
1673 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
1674 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
1675 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
1676 trk_px.clear(); trk_py.clear(); trk_pz.clear();
1677 trk_px.reserve(V0_helper.nRefTrks());
1678 trk_py.reserve(V0_helper.nRefTrks());
1679 trk_pz.reserve(V0_helper.nRefTrks());
1680 for(auto&& vec3 : V0_helper.refTrks()) {
1681 trk_px.push_back( vec3.Px() );
1682 trk_py.push_back( vec3.Py() );
1683 trk_pz.push_back( vec3.Pz() );
1684 }
1685 trk_pxDeco(*cascadeVertices[0]) = trk_px;
1686 trk_pyDeco(*cascadeVertices[0]) = trk_py;
1687 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
1688
1689 result.push_back( std::make_pair(fit_result.release(),nullptr) );
1690 }
1691 }
1692 } // loop over trackContainer
1693 } // for m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo<=0
1694 else if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0 && m_extraTrk3MassHypo<=0) {
1695 std::vector<const xAOD::TrackParticle*> tracksPlus;
1696 std::vector<const xAOD::TrackParticle*> tracksMinus;
1697 for(const xAOD::TrackParticle* tpExtra : *trackContainer) {
1698 if( tpExtra->pt() < std::fmin(m_extraTrk1MinPt,m_extraTrk2MinPt) ) continue;
1699 if( !m_trkSelector->decision(*tpExtra, nullptr) ) continue;
1700 // Check identical tracks in input
1701 if(std::find(tracksJX.cbegin(),tracksJX.cend(),tpExtra) != tracksJX.cend()) continue;
1702 if(std::find(tracksV0.cbegin(),tracksV0.cend(),tpExtra) != tracksV0.cend()) continue;
1703 if(tpExtra->charge()>0) {
1704 tracksPlus.push_back(tpExtra);
1705 }
1706 else {
1707 tracksMinus.push_back(tpExtra);
1708 }
1709 }
1710
1711 MesonCandidateVector D0Candidates(m_maxMesonCandidates, m_MesonPtOrdering);
1712 TLorentzVector p4_ExtraTrk1, p4_ExtraTrk2;
1713 for(const xAOD::TrackParticle* tp1 : tracksPlus) {
1714 for(const xAOD::TrackParticle* tp2 : tracksMinus) {
1715 if((tp1->pt()>m_extraTrk1MinPt && tp2->pt()>m_extraTrk2MinPt) ||
1716 (tp1->pt()>m_extraTrk2MinPt && tp2->pt()>m_extraTrk1MinPt)) {
1717 p4_ExtraTrk1.SetPtEtaPhiM(tp1->pt(), tp1->eta(), tp1->phi(), m_extraTrk1MassHypo);
1718 p4_ExtraTrk2.SetPtEtaPhiM(tp2->pt(), tp2->eta(), tp2->phi(), m_extraTrk2MassHypo);
1719 double main_mass = (p4_moth+p4_ExtraTrk1+p4_ExtraTrk2).M();
1720 if(m_useImprovedMass) {
1721 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_v0).M() + m_massJpsi;
1722 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_v0).M() + m_massJX;
1723 if((V0==LAMBDA || V0==LAMBDABAR) && m_massLd>0) main_mass += - p4_v0.M() + m_massLd;
1724 else if(V0==KS && m_massKs>0) main_mass += - p4_v0.M() + m_massKs;
1725 if(m_massD0>0) main_mass += - (p4_ExtraTrk1+p4_ExtraTrk2).M() + m_massD0;
1726 }
1727 if(main_mass < m_MassLower || main_mass > m_MassUpper) continue;
1728 auto D0 = getD0Candidate(ctx,JXvtx,tp1,tp2);
1729 if(D0.extraTrack1) D0Candidates.push_back(D0);
1730 }
1731 }
1732 }
1733
1734 std::vector<double> massesExtra{m_extraTrk1MassHypo,m_extraTrk2MassHypo};
1735
1736 for(auto&& D0 : D0Candidates.vector()) {
1737 std::vector<const xAOD::TrackParticle*> tracksExtra{D0.extraTrack1,D0.extraTrack2};
1738
1739 // Apply the user's settings to the fitter
1740 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
1741 // Robustness: http://cdsweb.cern.ch/record/685551
1742 int robustness = 0;
1743 m_iVertexFitter->setRobustness(robustness, *state);
1744 // Build up the topology
1745 // Vertex list
1746 std::vector<Trk::VertexID> vrtList;
1747 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
1748 // V0 vertex
1749 Trk::VertexID vID1;
1750 if (m_constrV0) {
1751 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,V0==KS?m_massKs:m_massLd);
1752 } else {
1753 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1754 }
1755 vrtList.push_back(vID1);
1756 // D0 vertex
1757 Trk::VertexID vID2;
1758 if (m_constrD0) {
1759 vID2 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state,m_massD0);
1760 } else {
1761 vID2 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state);
1762 }
1763 vrtList.push_back(vID2);
1764 // Mother vertex includes one subvertex (V0) and JX tracks + extra track
1765 Trk::VertexID vID3;
1766 if(m_constrMainV) {
1767 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state,m_massMainV);
1768 } else {
1769 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
1770 }
1771 if (m_constrJX && m_jxDaug_num>2) {
1772 std::vector<Trk::VertexID> cnstV;
1773 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
1774 ATH_MSG_WARNING("addMassConstraint for JX failed");
1775 }
1776 }
1777 if (m_constrJpsi) {
1778 std::vector<Trk::VertexID> cnstV;
1779 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
1780 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
1781 }
1782 }
1783 if (m_constrX && m_jxDaug_num==4) {
1784 std::vector<Trk::VertexID> cnstV;
1785 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksX,cnstV,*state,m_massX).isSuccess() ) {
1786 ATH_MSG_WARNING("addMassConstraint for X failed");
1787 }
1788 }
1789 // Do the work
1790 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) );
1791
1792 if (fit_result) {
1793 for(auto& v : fit_result->vertices()) {
1794 if(v->nTrackParticles()==0) {
1795 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
1796 v->setTrackParticleLinks(nullLinkVector);
1797 }
1798 }
1799 // reset links to original tracks
1800 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
1801
1802 // necessary to prevent memory leak
1803 fit_result->setSVOwnership(true);
1804
1805 // Chi2/DOF cut
1806 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
1807 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
1808 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
1809 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
1810 size_t iMoth = cascadeVertices.size()-1;
1811 double lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[iMoth]);
1812 double lxy_SV2 = m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[iMoth]);
1813 if(chi2CutPassed && lxy_SV1>m_lxyV0_cut && lxy_SV2>m_lxyD0_cut) {
1814 chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
1815 ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
1816 if(V0==LAMBDA) {
1817 type_V1_decor(*cascadeVertices[0]) = "Lambda";
1818 }
1819 else if(V0==LAMBDABAR) {
1820 type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
1821 }
1822 else if(V0==KS) {
1823 type_V1_decor(*cascadeVertices[0]) = "Ks";
1824 }
1825 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
1826 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
1827 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
1828 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
1829 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
1830 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
1831 trk_px.clear(); trk_py.clear(); trk_pz.clear();
1832 trk_px.reserve(V0_helper.nRefTrks());
1833 trk_py.reserve(V0_helper.nRefTrks());
1834 trk_pz.reserve(V0_helper.nRefTrks());
1835 for(auto&& vec3 : V0_helper.refTrks()) {
1836 trk_px.push_back( vec3.Px() );
1837 trk_py.push_back( vec3.Py() );
1838 trk_pz.push_back( vec3.Pz() );
1839 }
1840 trk_pxDeco(*cascadeVertices[0]) = trk_px;
1841 trk_pyDeco(*cascadeVertices[0]) = trk_py;
1842 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
1843
1844 result.push_back( std::make_pair(fit_result.release(),nullptr) );
1845 }
1846 }
1847 } // loop over D0Candidates
1848 } // for m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0 && m_extraTrk3MassHypo<=0
1849 else if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0 && m_extraTrk3MassHypo>0) {
1850 std::vector<const xAOD::TrackParticle*> tracksPlus;
1851 std::vector<const xAOD::TrackParticle*> tracksMinus;
1852 double minTrkPt = std::fmin(std::fmin(m_extraTrk1MinPt,m_extraTrk2MinPt),m_extraTrk3MinPt);
1853 for(const xAOD::TrackParticle* tpExtra : *trackContainer) {
1854 if( tpExtra->pt() < minTrkPt ) continue;
1855 if( !m_trkSelector->decision(*tpExtra, nullptr) ) continue;
1856 // Check identical tracks in input
1857 if(std::find(tracksJX.cbegin(),tracksJX.cend(),tpExtra) != tracksJX.cend()) continue;
1858 if(std::find(tracksV0.cbegin(),tracksV0.cend(),tpExtra) != tracksV0.cend()) continue;
1859 if(tpExtra->charge()>0) {
1860 tracksPlus.push_back(tpExtra);
1861 }
1862 else {
1863 tracksMinus.push_back(tpExtra);
1864 }
1865 }
1866
1867 MesonCandidateVector DpmCandidates(m_maxMesonCandidates, m_MesonPtOrdering);
1868 TLorentzVector p4_ExtraTrk1, p4_ExtraTrk2, p4_ExtraTrk3;
1869 // +,-,- combination (D- -> K+ pi- pi-)
1870 for(const xAOD::TrackParticle* tp1 : tracksPlus) {
1871 if( tp1->pt() < m_extraTrk1MinPt ) continue;
1872 for(auto tp2Itr=tracksMinus.cbegin(); tp2Itr!=tracksMinus.cend(); ++tp2Itr) {
1873 const xAOD::TrackParticle* tp2 = *tp2Itr;
1874 for(auto tp3Itr=tp2Itr+1; tp3Itr!=tracksMinus.cend(); ++tp3Itr) {
1875 const xAOD::TrackParticle* tp3 = *tp3Itr;
1876 if((tp2->pt()>m_extraTrk2MinPt && tp3->pt()>m_extraTrk3MinPt) ||
1877 (tp2->pt()>m_extraTrk3MinPt && tp3->pt()>m_extraTrk2MinPt)) {
1878 p4_ExtraTrk1.SetPtEtaPhiM(tp1->pt(), tp1->eta(), tp1->phi(), m_extraTrk1MassHypo);
1879 p4_ExtraTrk2.SetPtEtaPhiM(tp2->pt(), tp2->eta(), tp2->phi(), m_extraTrk2MassHypo);
1880 p4_ExtraTrk3.SetPtEtaPhiM(tp3->pt(), tp3->eta(), tp3->phi(), m_extraTrk3MassHypo);
1881 double main_mass = (p4_moth+p4_ExtraTrk1+p4_ExtraTrk2+p4_ExtraTrk3).M();
1882 if(m_useImprovedMass) {
1883 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_v0).M() + m_massJpsi;
1884 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_v0).M() + m_massJX;
1885 if((V0==LAMBDA || V0==LAMBDABAR) && m_massLd>0) main_mass += - p4_v0.M() + m_massLd;
1886 else if(V0==KS && m_massKs>0) main_mass += - p4_v0.M() + m_massKs;
1887 if(m_massDpm>0) main_mass += - (p4_ExtraTrk1+p4_ExtraTrk2+p4_ExtraTrk3).M() + m_massDpm;
1888 }
1889 if(main_mass < m_MassLower || main_mass > m_MassUpper) continue;
1890 auto Dpm = getDpmCandidate(ctx,JXvtx,tp1,tp2,tp3);
1891 if(Dpm.extraTrack1) DpmCandidates.push_back(Dpm);
1892 }
1893 }
1894 }
1895 }
1896 // -,+,+ combination (D+ -> K- pi+ pi+)
1897 for(const xAOD::TrackParticle* tp1 : tracksMinus) {
1898 if( tp1->pt() < m_extraTrk1MinPt ) continue;
1899 for(auto tp2Itr=tracksPlus.cbegin(); tp2Itr!=tracksPlus.cend(); ++tp2Itr) {
1900 const xAOD::TrackParticle* tp2 = *tp2Itr;
1901 for(auto tp3Itr=tp2Itr+1; tp3Itr!=tracksPlus.cend(); ++tp3Itr) {
1902 const xAOD::TrackParticle* tp3 = *tp3Itr;
1903 if((tp2->pt()>m_extraTrk2MinPt && tp3->pt()>m_extraTrk3MinPt) ||
1904 (tp2->pt()>m_extraTrk3MinPt && tp3->pt()>m_extraTrk2MinPt)) {
1905 p4_ExtraTrk1.SetPtEtaPhiM(tp1->pt(), tp1->eta(), tp1->phi(), m_extraTrk1MassHypo);
1906 p4_ExtraTrk2.SetPtEtaPhiM(tp2->pt(), tp2->eta(), tp2->phi(), m_extraTrk2MassHypo);
1907 p4_ExtraTrk3.SetPtEtaPhiM(tp3->pt(), tp3->eta(), tp3->phi(), m_extraTrk3MassHypo);
1908 double main_mass = (p4_moth+p4_ExtraTrk1+p4_ExtraTrk2+p4_ExtraTrk3).M();
1909 if(m_useImprovedMass) {
1910 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_v0).M() + m_massJpsi;
1911 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_v0).M() + m_massJX;
1912 if((V0==LAMBDA || V0==LAMBDABAR) && m_massLd>0) main_mass += - p4_v0.M() + m_massLd;
1913 else if(V0==KS && m_massKs>0) main_mass += - p4_v0.M() + m_massKs;
1914 if(m_massDpm>0) main_mass += - (p4_ExtraTrk1+p4_ExtraTrk2+p4_ExtraTrk3).M() + m_massDpm;
1915 }
1916 if(main_mass < m_MassLower || main_mass > m_MassUpper) continue;
1917 auto Dpm = getDpmCandidate(ctx,JXvtx,tp1,tp2,tp3);
1918 if(Dpm.extraTrack1) DpmCandidates.push_back(Dpm);
1919 }
1920 }
1921 }
1922 }
1923
1924 std::vector<double> massesExtra{m_extraTrk1MassHypo,m_extraTrk2MassHypo,m_extraTrk3MassHypo};
1925
1926 for(auto&& Dpm : DpmCandidates.vector()) {
1927 std::vector<const xAOD::TrackParticle*> tracksExtra{Dpm.extraTrack1,Dpm.extraTrack2,Dpm.extraTrack3};
1928
1929 // Apply the user's settings to the fitter
1930 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
1931 // Robustness: http://cdsweb.cern.ch/record/685551
1932 int robustness = 0;
1933 m_iVertexFitter->setRobustness(robustness, *state);
1934 // Build up the topology
1935 // Vertex list
1936 std::vector<Trk::VertexID> vrtList;
1937 std::vector<Trk::VertexID> vrtList2;
1938 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
1939 if(m_JXV0SubVtx) {
1940 // V0 vertex
1941 Trk::VertexID vID1;
1942 if (m_constrV0) {
1943 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,V0==KS?m_massKs:m_massLd);
1944 } else {
1945 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1946 }
1947 vrtList.push_back(vID1);
1948 // JX+V0 vertex
1949 Trk::VertexID vID2;
1950 if(m_constrJXV0) {
1951 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state,m_massJXV0);
1952 } else {
1953 vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
1954 }
1955 vrtList2.push_back(vID2);
1956 Trk::VertexID vID3;
1957 if (m_constrDpm) {
1958 vID3 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state,m_massDpm);
1959 } else {
1960 vID3 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state);
1961 }
1962 vrtList2.push_back(vID3);
1963 // Mother vertex includes two subvertices: V0, Dpm and JX tracks
1964 std::vector<const xAOD::TrackParticle*> tp;
1965 std::vector<double> tp_masses;
1966 if(m_constrMainV) {
1967 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2,*state,m_massMainV);
1968 } else {
1969 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2,*state);
1970 }
1971 if (m_constrJX && m_jxDaug_num>2) {
1972 std::vector<Trk::VertexID> cnstV;
1973 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
1974 ATH_MSG_WARNING("addMassConstraint for JX failed");
1975 }
1976 }
1977 if (m_constrJpsi) {
1978 std::vector<Trk::VertexID> cnstV;
1979 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
1980 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
1981 }
1982 }
1983 if (m_constrX && m_jxDaug_num==4) {
1984 std::vector<Trk::VertexID> cnstV;
1985 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksX,cnstV,*state,m_massX).isSuccess() ) {
1986 ATH_MSG_WARNING("addMassConstraint for X failed");
1987 }
1988 }
1989 }
1990 else { // m_JXV0SubVtx==false
1991 // V0 vertex
1992 Trk::VertexID vID1;
1993 if (m_constrV0) {
1994 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,V0==KS?m_massKs:m_massLd);
1995 } else {
1996 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1997 }
1998 vrtList.push_back(vID1);
1999 // Dpm vertex
2000 Trk::VertexID vID2;
2001 if (m_constrDpm) {
2002 vID2 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state,m_massDpm);
2003 } else {
2004 vID2 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state);
2005 }
2006 vrtList.push_back(vID2);
2007 // Mother vertex includes two subvertices: V0, Dpm and JX tracks
2008 Trk::VertexID vID3;
2009 if(m_constrMainV) {
2010 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state,m_massMainV);
2011 } else {
2012 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
2013 }
2014 if (m_constrJX && m_jxDaug_num>2) {
2015 std::vector<Trk::VertexID> cnstV;
2016 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
2017 ATH_MSG_WARNING("addMassConstraint for JX failed");
2018 }
2019 }
2020 if (m_constrJpsi) {
2021 std::vector<Trk::VertexID> cnstV;
2022 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
2023 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
2024 }
2025 }
2026 if (m_constrX && m_jxDaug_num==4) {
2027 std::vector<Trk::VertexID> cnstV;
2028 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksX,cnstV,*state,m_massX).isSuccess() ) {
2029 ATH_MSG_WARNING("addMassConstraint for X failed");
2030 }
2031 }
2032 }
2033 // Do the work
2034 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) );
2035
2036 if (fit_result) {
2037 for(auto& v : fit_result->vertices()) {
2038 if(v->nTrackParticles()==0) {
2039 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
2040 v->setTrackParticleLinks(nullLinkVector);
2041 }
2042 }
2043 // reset links to original tracks
2044 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
2045
2046 // necessary to prevent memory leak
2047 fit_result->setSVOwnership(true);
2048
2049 // Chi2/DOF cut
2050 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
2051 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
2052 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
2053 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
2054 size_t iMoth = cascadeVertices.size()-1;
2055 double lxy_SV1(0), lxy_SV2(0);
2056 if(m_JXV0SubVtx) {
2057 lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
2058 lxy_SV2 = m_CascadeTools->lxy(moms[2],cascadeVertices[2],cascadeVertices[iMoth]);
2059 }
2060 else {
2061 lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[iMoth]);
2062 lxy_SV2 = m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[iMoth]);
2063 }
2064 if(chi2CutPassed && lxy_SV1>m_lxyV0_cut && lxy_SV2>m_lxyDpm_cut) {
2065 chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
2066 ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
2067 if(V0==LAMBDA) {
2068 type_V1_decor(*cascadeVertices[0]) = "Lambda";
2069 }
2070 else if(V0==LAMBDABAR) {
2071 type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
2072 }
2073 else if(V0==KS) {
2074 type_V1_decor(*cascadeVertices[0]) = "Ks";
2075 }
2076 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
2077 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
2078 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
2079 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
2080 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
2081 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
2082 trk_px.clear(); trk_py.clear(); trk_pz.clear();
2083 trk_px.reserve(V0_helper.nRefTrks());
2084 trk_py.reserve(V0_helper.nRefTrks());
2085 trk_pz.reserve(V0_helper.nRefTrks());
2086 for(auto&& vec3 : V0_helper.refTrks()) {
2087 trk_px.push_back( vec3.Px() );
2088 trk_py.push_back( vec3.Py() );
2089 trk_pz.push_back( vec3.Pz() );
2090 }
2091 trk_pxDeco(*cascadeVertices[0]) = trk_px;
2092 trk_pyDeco(*cascadeVertices[0]) = trk_py;
2093 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
2094
2095 // refit with main vertex mass constraint if required
2096 if(!m_constrMainV && m_doPostMainVContrFit) {
2097 TLorentzVector totalMom;
2098 for(size_t it=0; it<moms[iMoth].size(); it++) totalMom += moms[iMoth][it];
2099 double mainV_mass = totalMom.M();
2100 if(mainV_mass>m_PostMassLower && mainV_mass<m_PostMassUpper) {
2101 std::unique_ptr<Trk::IVKalState> state_mvc = m_iVertexFitter->makeState(ctx);
2102 int robustness_mvc = 0;
2103 m_iVertexFitter->setRobustness(robustness_mvc, *state_mvc);
2104 std::vector<Trk::VertexID> vrtList_mvc;
2105 std::vector<Trk::VertexID> vrtList2_mvc;
2106 if(m_JXV0SubVtx) {
2107 Trk::VertexID vID1_mvc; // V0 vertex
2108 if (m_constrV0) {
2109 vID1_mvc = m_iVertexFitter->startVertex(tracksV0,massesV0,*state_mvc,V0==KS?m_massKs:m_massLd);
2110 } else {
2111 vID1_mvc = m_iVertexFitter->startVertex(tracksV0,massesV0,*state_mvc);
2112 }
2113 vrtList_mvc.push_back(vID1_mvc);
2114 Trk::VertexID vID2_mvc; // JX+V0 vertex
2115 if(m_constrJXV0) {
2116 vID2_mvc = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList_mvc,*state_mvc,m_massJXV0);
2117 } else {
2118 vID2_mvc = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList_mvc,*state_mvc);
2119 }
2120 vrtList2_mvc.push_back(vID2_mvc);
2121 Trk::VertexID vID3_mvc; // Dpm vertex
2122 if (m_constrDpm) {
2123 vID3_mvc = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state_mvc,m_massDpm);
2124 } else {
2125 vID3_mvc = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state_mvc);
2126 }
2127 vrtList2_mvc.push_back(vID3_mvc);
2128 std::vector<const xAOD::TrackParticle*> tp;
2129 std::vector<double> tp_masses;
2130 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2_mvc,*state_mvc,m_massMainV); // Mother vertex
2131 if (m_constrJX && m_jxDaug_num>2) {
2132 std::vector<Trk::VertexID> cnstV_mvc;
2133 if ( !m_iVertexFitter->addMassConstraint(vID2_mvc,tracksJX,cnstV_mvc,*state_mvc,m_massJX).isSuccess() ) {
2134 ATH_MSG_WARNING("addMassConstraint for JX failed");
2135 }
2136 }
2137 if (m_constrJpsi) {
2138 std::vector<Trk::VertexID> cnstV_mvc;
2139 if ( !m_iVertexFitter->addMassConstraint(vID2_mvc,tracksJpsi,cnstV_mvc,*state_mvc,m_massJpsi).isSuccess() ) {
2140 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
2141 }
2142 }
2143 if (m_constrX && m_jxDaug_num==4) {
2144 std::vector<Trk::VertexID> cnstV_mvc;
2145 if ( !m_iVertexFitter->addMassConstraint(vID2_mvc,tracksX,cnstV_mvc,*state_mvc,m_massX).isSuccess() ) {
2146 ATH_MSG_WARNING("addMassConstraint for X failed");
2147 }
2148 }
2149 }
2150 else {
2151 Trk::VertexID vID1_mvc; // V0 vertex
2152 if (m_constrV0) {
2153 vID1_mvc = m_iVertexFitter->startVertex(tracksV0,massesV0,*state_mvc,V0==KS?m_massKs:m_massLd);
2154 } else {
2155 vID1_mvc = m_iVertexFitter->startVertex(tracksV0,massesV0,*state_mvc);
2156 }
2157 vrtList_mvc.push_back(vID1_mvc);
2158 Trk::VertexID vID2_mvc; // Dpm vertex
2159 if (m_constrDpm) {
2160 vID2_mvc = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state_mvc,m_massDpm);
2161 } else {
2162 vID2_mvc = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state_mvc);
2163 }
2164 vrtList_mvc.push_back(vID2_mvc);
2165 Trk::VertexID vID3_mvc = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList_mvc,*state_mvc,m_massMainV); // Mother vertex
2166 if (m_constrJX && m_jxDaug_num>2) {
2167 std::vector<Trk::VertexID> cnstV_mvc;
2168 if ( !m_iVertexFitter->addMassConstraint(vID3_mvc,tracksJX,cnstV_mvc,*state_mvc,m_massJX).isSuccess() ) {
2169 ATH_MSG_WARNING("addMassConstraint for JX failed");
2170 }
2171 }
2172 if (m_constrJpsi) {
2173 std::vector<Trk::VertexID> cnstV_mvc;
2174 if ( !m_iVertexFitter->addMassConstraint(vID3_mvc,tracksJpsi,cnstV_mvc,*state_mvc,m_massJpsi).isSuccess() ) {
2175 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
2176 }
2177 }
2178 if (m_constrX && m_jxDaug_num==4) {
2179 std::vector<Trk::VertexID> cnstV_mvc;
2180 if ( !m_iVertexFitter->addMassConstraint(vID3_mvc,tracksX,cnstV_mvc,*state_mvc,m_massX).isSuccess() ) {
2181 ATH_MSG_WARNING("addMassConstraint for X failed");
2182 }
2183 }
2184 }
2185 // Do the work
2186 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) );
2187
2188 if (fit_result_mvc) {
2189 for(auto& v : fit_result_mvc->vertices()) {
2190 if(v->nTrackParticles()==0) {
2191 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
2192 v->setTrackParticleLinks(nullLinkVector);
2193 }
2194 }
2195 BPhysPVCascadeTools::PrepareVertexLinks(fit_result_mvc.get(), trackCols);
2196 fit_result_mvc->setSVOwnership(true);
2197 chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
2198 ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
2199 if(V0==LAMBDA) {
2200 type_V1_decor(*cascadeVertices[0]) = "Lambda";
2201 }
2202 else if(V0==LAMBDABAR) {
2203 type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
2204 }
2205 else if(V0==KS) {
2206 type_V1_decor(*cascadeVertices[0]) = "Ks";
2207 }
2208 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
2209 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
2210 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
2211 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
2212 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
2213 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
2214 trk_px.clear(); trk_py.clear(); trk_pz.clear();
2215 trk_px.reserve(V0_helper.nRefTrks());
2216 trk_py.reserve(V0_helper.nRefTrks());
2217 trk_pz.reserve(V0_helper.nRefTrks());
2218 for(auto&& vec3 : V0_helper.refTrks()) {
2219 trk_px.push_back( vec3.Px() );
2220 trk_py.push_back( vec3.Py() );
2221 trk_pz.push_back( vec3.Pz() );
2222 }
2223 trk_pxDeco(*cascadeVertices[0]) = trk_px;
2224 trk_pyDeco(*cascadeVertices[0]) = trk_py;
2225 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
2226
2227 result.push_back( std::make_pair(fit_result.release(),fit_result_mvc.release()) );
2228 }
2229 else result.push_back( std::make_pair(fit_result.release(),nullptr) );
2230 }
2231 else result.push_back( std::make_pair(fit_result.release(),nullptr) );
2232 }
2233 else result.push_back( std::make_pair(fit_result.release(),nullptr) );
2234 }
2235 }
2236 } // loop over DpmCandidates
2237 } // for m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0 && m_extraTrk3MassHypo>0
2238
2239 if(pv_xAOD) {
2240 for(auto cascade_info_pair : result) {
2241 if(cascade_info_pair.first && cascade_info_pair.first->getParticleMoms().size()>0) {
2242 size_t index = cascade_info_pair.first->getParticleMoms().size() - 1;
2243 const std::vector<TLorentzVector> &mom = cascade_info_pair.first->getParticleMoms()[index];
2244 const Amg::MatrixX &cov = cascade_info_pair.first->getCovariance()[index];
2245 const xAOD::Vertex* mainVertex = cascade_info_pair.first->vertices()[index];
2246 xAOD::BPhysHypoHelper vtx(m_hypoName, mainVertex);
2247 bool isInDefaultPVCont = false;
2248 for(const xAOD::Vertex* pvVtx : *defaultPVContainer) {
2249 if(pv_xAOD == pvVtx) { isInDefaultPVCont = true; break; }
2250 }
2251 if(isInDefaultPVCont) vtx.setPv( pv_xAOD, defaultPVContainer, pvtype );
2252 else vtx.setPv( pv_xAOD, pvContainer, pvtype );
2253 if(origPv_xAOD) vtx.setOrigPv( origPv_xAOD, defaultPVContainer, pvtype );
2254 vtx.setLxy ( m_CascadeTools->lxy (mom, vtx.vtx(), pv_xAOD), pvtype );
2255 vtx.setLxyErr ( m_CascadeTools->lxyError (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2256 vtx.setA0 ( m_CascadeTools->a0 (mom, vtx.vtx(), pv_xAOD), pvtype );
2257 vtx.setA0Err ( m_CascadeTools->a0Error (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2258 vtx.setA0xy ( m_CascadeTools->a0xy (mom, vtx.vtx(), pv_xAOD), pvtype );
2259 vtx.setA0xyErr( m_CascadeTools->a0xyError(mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2260 vtx.setZ0 ( m_CascadeTools->a0z (mom, vtx.vtx(), pv_xAOD), pvtype );
2261 vtx.setZ0Err ( m_CascadeTools->a0zError (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2262 vtx.setRefitPVStatus( 0, pvtype );
2263 // Proper decay times
2264 vtx.setTau( m_CascadeTools->tau(mom, vtx.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
2265 vtx.setTauErr( m_CascadeTools->tauError(mom, cov, vtx.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
2266 vtx.setTau( m_CascadeTools->tau(mom, vtx.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
2267 vtx.setTauErr( m_CascadeTools->tauError(mom, cov, vtx.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
2268
2269 if(cascade_info_pair.second && cascade_info_pair.second->getParticleMoms().size()>0) {
2270 index = cascade_info_pair.second->getParticleMoms().size() - 1;
2271 const std::vector<TLorentzVector> &mom_mvc = cascade_info_pair.second->getParticleMoms()[index];
2272 const Amg::MatrixX &cov_mvc = cascade_info_pair.second->getCovariance()[index];
2273 const xAOD::Vertex* mainVertex_mvc = cascade_info_pair.second->vertices()[index];
2274 xAOD::BPhysHypoHelper vtx_mvc(m_hypoName, mainVertex_mvc);
2275 if(isInDefaultPVCont) vtx_mvc.setPv( pv_xAOD, defaultPVContainer, pvtype );
2276 else vtx_mvc.setPv( pv_xAOD, pvContainer, pvtype );
2277 if(origPv_xAOD) vtx.setOrigPv( origPv_xAOD, defaultPVContainer, pvtype );
2278 vtx_mvc.setLxy ( m_CascadeTools->lxy (mom_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2279 vtx_mvc.setLxyErr ( m_CascadeTools->lxyError (mom_mvc, cov_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2280 vtx_mvc.setA0 ( m_CascadeTools->a0 (mom_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2281 vtx_mvc.setA0Err ( m_CascadeTools->a0Error (mom_mvc, cov_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2282 vtx_mvc.setA0xy ( m_CascadeTools->a0xy (mom_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2283 vtx_mvc.setA0xyErr( m_CascadeTools->a0xyError(mom_mvc, cov_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2284 vtx_mvc.setZ0 ( m_CascadeTools->a0z (mom_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2285 vtx_mvc.setZ0Err ( m_CascadeTools->a0zError (mom_mvc, cov_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype );
2286 vtx_mvc.setRefitPVStatus( 0, pvtype );
2287 // Proper decay times
2288 vtx_mvc.setTau( m_CascadeTools->tau(mom_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
2289 vtx_mvc.setTauErr( m_CascadeTools->tauError(mom_mvc, cov_mvc, vtx_mvc.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
2290 vtx_mvc.setTau( m_CascadeTools->tau(mom_mvc, vtx_mvc.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
2291 vtx_mvc.setTauErr( m_CascadeTools->tauError(mom_mvc, cov_mvc, vtx_mvc.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
2292 }
2293 }
2294 }
2295 }
2296
2297 return result;
2298 }
2299
2300 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 {
2301 std::vector<std::pair<Trk::VxCascadeInfo*,Trk::VxCascadeInfo*> > result;
2302
2303 std::vector<const xAOD::TrackParticle*> tracksJX;
2304 tracksJX.reserve(JXvtx->nTrackParticles());
2305 for(size_t i=0; i<JXvtx->nTrackParticles(); i++) tracksJX.push_back(JXvtx->trackParticle(i));
2306 if (tracksJX.size() != massesJX.size()) {
2307 ATH_MSG_ERROR("Problems with JX input: number of tracks or track mass inputs is not correct!");
2308 return result;
2309 }
2310 // Check identical tracks in input
2311 if(std::find(tracksJX.cbegin(), tracksJX.cend(), disVtx.V0vtx->trackParticle(0)) != tracksJX.cend()) return result;
2312 if(std::find(tracksJX.cbegin(), tracksJX.cend(), disVtx.V0vtx->trackParticle(1)) != tracksJX.cend()) return result;
2313 std::vector<const xAOD::TrackParticle*> tracksV0;
2314 tracksV0.reserve(disVtx.V0vtx->nTrackParticles());
2315 for(size_t j=0; j<disVtx.V0vtx->nTrackParticles(); j++) tracksV0.push_back(disVtx.V0vtx->trackParticle(j));
2316
2317 if(std::find(tracksJX.cbegin(), tracksJX.cend(), disVtx.track) != tracksJX.cend()) return result;
2318 std::vector<const xAOD::TrackParticle*> tracks3{disVtx.track};
2319 std::vector<double> massesDis3{m_disVDaug3MassHypo};
2320
2321 std::vector<const xAOD::TrackParticle*> tracksJpsi{tracksJX[0], tracksJX[1]};
2322 std::vector<const xAOD::TrackParticle*> tracksX;
2323 if(m_jxDaug_num>=3) tracksX.push_back(tracksJX[2]);
2324 if(m_jxDaug_num==4) tracksX.push_back(tracksJX[3]);
2325
2326 std::vector<double> massesV0;
2327 if(disVtx.V0type==LAMBDA) {
2328 massesV0 = m_massesV0_ppi;
2329 }
2330 else if(disVtx.V0type==LAMBDABAR) {
2331 massesV0 = m_massesV0_pip;
2332 }
2333 else if(disVtx.V0type==KS) {
2334 massesV0 = m_massesV0_pipi;
2335 }
2336
2337 std::vector<double> massesDisV = massesV0; massesDisV.push_back(m_disVDaug3MassHypo);
2338
2339 TLorentzVector p4_moth, p4_disV, tmp;
2340 for(size_t it=0; it<JXvtx->nTrackParticles(); it++) {
2341 tmp.SetPtEtaPhiM(JXvtx->trackParticle(it)->pt(), JXvtx->trackParticle(it)->eta(), JXvtx->trackParticle(it)->phi(), massesJX[it]);
2342 p4_moth += tmp;
2343 }
2344 p4_disV += disVtx.p4_V0track1; p4_disV += disVtx.p4_V0track2; p4_disV += disVtx.p4_disVtrack;
2345 p4_moth += p4_disV;
2346
2347 xAOD::BPhysHelper::pv_type pvtype = xAOD::BPhysHelper::PV_MIN_A0;
2348 if(m_cascadeFitWithPV==1) pvtype = xAOD::BPhysHelper::PV_MAX_SUM_PT2;
2349 else if(m_cascadeFitWithPV==2) pvtype = xAOD::BPhysHelper::PV_MIN_A0;
2350 else if(m_cascadeFitWithPV==3) pvtype = xAOD::BPhysHelper::PV_MIN_Z0;
2351 else if(m_cascadeFitWithPV==4) pvtype = xAOD::BPhysHelper::PV_MIN_Z0_BA;
2352 xAOD::BPhysHelper JX_helper(JXvtx);
2353 const xAOD::Vertex* origPv_xAOD = (m_cascadeFitWithPV>=1 && m_cascadeFitWithPV<=4) ? JX_helper.origPv(pvtype) : nullptr;
2354 const xAOD::Vertex* pv_xAOD = (m_cascadeFitWithPV>=1 && m_cascadeFitWithPV<=4) ? JX_helper.pv(pvtype) : nullptr;
2355 std::unique_ptr<Trk::RecVertex> pv_AOD;
2356 if(pv_xAOD) pv_AOD = std::make_unique<Trk::RecVertex>(pv_xAOD->position(),pv_xAOD->covariancePosition(),pv_xAOD->numberDoF(),pv_xAOD->chiSquared());
2357
2358 SG::Decorator<float> chi2_V0_decor("ChiSquared_V0");
2359 SG::Decorator<int> ndof_V0_decor("nDoF_V0");
2360 SG::Decorator<std::string> type_V0_decor("Type_V0");
2361
2362 SG::Accessor<int> mAcc_gfit("gamma_fit");
2363 SG::Accessor<float> mAcc_gmass("gamma_mass");
2364 SG::Accessor<float> mAcc_gmasserr("gamma_massError");
2365 SG::Accessor<float> mAcc_gchisq("gamma_chisq");
2366 SG::Accessor<int> mAcc_gndof("gamma_ndof");
2367 SG::Accessor<float> mAcc_gprob("gamma_probability");
2368
2369 SG::Decorator<int> mDec_gfit("gamma_fit");
2370 SG::Decorator<float> mDec_gmass("gamma_mass");
2371 SG::Decorator<float> mDec_gmasserr("gamma_massError");
2372 SG::Decorator<float> mDec_gchisq("gamma_chisq");
2373 SG::Decorator<int> mDec_gndof("gamma_ndof");
2374 SG::Decorator<float> mDec_gprob("gamma_probability");
2375 SG::Decorator< std::vector<float> > trk_pxDeco("TrackPx_V0nc");
2376 SG::Decorator< std::vector<float> > trk_pyDeco("TrackPy_V0nc");
2377 SG::Decorator< std::vector<float> > trk_pzDeco("TrackPz_V0nc");
2378 SG::Decorator<float> trk_px_deco("TrackPx_DisVnc");
2379 SG::Decorator<float> trk_py_deco("TrackPy_DisVnc");
2380 SG::Decorator<float> trk_pz_deco("TrackPz_DisVnc");
2381
2382 std::vector<float> trk_px;
2383 std::vector<float> trk_py;
2384 std::vector<float> trk_pz;
2385
2386 if(m_extraTrk1MassHypo<=0) {
2387 double main_mass = p4_moth.M();
2388 if(m_useImprovedMass) {
2389 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_disV).M() + m_massJpsi;
2390 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_disV).M() + m_massJX;
2391 if(m_massDisV>0) main_mass += - p4_disV.M() + m_massDisV;
2392 }
2393 if (main_mass < m_MassLower || main_mass > m_MassUpper) return result;
2394
2395 // Apply the user's settings to the fitter
2396 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
2397 // Robustness: http://cdsweb.cern.ch/record/685551
2398 int robustness = 0;
2399 m_iVertexFitter->setRobustness(robustness, *state);
2400 // Build up the topology
2401 // Vertex list
2402 std::vector<Trk::VertexID> vrtList;
2403 std::vector<Trk::VertexID> vrtList2;
2404 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
2405 // V0 vertex
2406 Trk::VertexID vID1;
2407 if (m_constrV0) {
2408 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,disVtx.V0type==KS?m_massKs:m_massLd);
2409 } else {
2410 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
2411 }
2412 vrtList.push_back(vID1);
2413 // Displaced vertex
2414 Trk::VertexID vID2;
2415 if (m_constrDisV) {
2416 vID2 = m_iVertexFitter->nextVertex(tracks3,massesDis3,vrtList,*state,m_massDisV);
2417 } else {
2418 vID2 = m_iVertexFitter->nextVertex(tracks3,massesDis3,vrtList,*state);
2419 }
2420 vrtList2.push_back(vID2);
2421 Trk::VertexID vID3;
2422 if(m_JXSubVtx) {
2423 // JX vertex
2424 if (m_constrJX && m_jxDaug_num>2) {
2425 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state,m_massJX);
2426 } else {
2427 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state);
2428 }
2429 vrtList2.push_back(vID3);
2430 // Mother vertex includes two subvertices: DisV and JX
2431 std::vector<const xAOD::TrackParticle*> tp;
2432 std::vector<double> tp_masses;
2433 if(m_constrMainV) {
2434 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2,*state,m_massMainV);
2435 } else {
2436 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2,*state);
2437 }
2438 }
2439 else { // m_JXSubVtx=false
2440 // Mother vertex includes just one subvertex (DisV) and JX tracks
2441 if(m_constrMainV) {
2442 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList2,*state,m_massMainV);
2443 } else {
2444 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList2,*state);
2445 }
2446 if (m_constrJX && m_jxDaug_num>2) {
2447 std::vector<Trk::VertexID> cnstV;
2448 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
2449 ATH_MSG_WARNING("addMassConstraint for JX failed");
2450 }
2451 }
2452 }
2453 if (m_constrJpsi) {
2454 std::vector<Trk::VertexID> cnstV;
2455 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
2456 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
2457 }
2458 }
2459 if (m_constrX && m_jxDaug_num==4) {
2460 std::vector<Trk::VertexID> cnstV;
2461 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksX,cnstV,*state,m_massX).isSuccess() ) {
2462 ATH_MSG_WARNING("addMassConstraint for X failed");
2463 }
2464 }
2465 // Do the work
2466 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) );
2467
2468 if (fit_result) {
2469 for(auto& v : fit_result->vertices()) {
2470 if(v->nTrackParticles()==0) {
2471 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
2472 v->setTrackParticleLinks(nullLinkVector);
2473 }
2474 }
2475 // reset links to original tracks
2476 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
2477
2478 // necessary to prevent memory leak
2479 fit_result->setSVOwnership(true);
2480
2481 // Chi2/DOF cut
2482 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
2483 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
2484
2485 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
2486 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
2487 size_t iMoth = cascadeVertices.size()-1;
2488 double lxy_SV1_sub = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
2489 double lxy_SV1 = m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[iMoth]);
2490
2491 if(chi2CutPassed && lxy_SV1>m_lxyDisV_cut && lxy_SV1_sub>m_lxyV0_cut) {
2492 chi2_V0_decor(*cascadeVertices[0]) = disVtx.V0vtx->chiSquared();
2493 ndof_V0_decor(*cascadeVertices[0]) = disVtx.V0vtx->numberDoF();
2494 if(disVtx.V0type==LAMBDA) {
2495 type_V0_decor(*cascadeVertices[0]) = "Lambda";
2496 }
2497 else if(disVtx.V0type==LAMBDABAR) {
2498 type_V0_decor(*cascadeVertices[0]) = "Lambdabar";
2499 }
2500 else if(disVtx.V0type==KS) {
2501 type_V0_decor(*cascadeVertices[0]) = "Ks";
2502 }
2503 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*disVtx.V0vtx) ? mAcc_gfit(*disVtx.V0vtx) : 0;
2504 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*disVtx.V0vtx) ? mAcc_gmass(*disVtx.V0vtx) : -1;
2505 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*disVtx.V0vtx) ? mAcc_gmasserr(*disVtx.V0vtx) : -1;
2506 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*disVtx.V0vtx) ? mAcc_gchisq(*disVtx.V0vtx) : 999999;
2507 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*disVtx.V0vtx) ? mAcc_gndof(*disVtx.V0vtx) : 0;
2508 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*disVtx.V0vtx) ? mAcc_gprob(*disVtx.V0vtx) : -1;
2509 trk_px.clear(); trk_py.clear(); trk_pz.clear();
2510 trk_px.push_back( disVtx.p4_V0track1.Px() ); trk_px.push_back( disVtx.p4_V0track2.Px() );
2511 trk_py.push_back( disVtx.p4_V0track1.Py() ); trk_py.push_back( disVtx.p4_V0track2.Py() );
2512 trk_pz.push_back( disVtx.p4_V0track1.Pz() ); trk_pz.push_back( disVtx.p4_V0track2.Pz() );
2513 trk_pxDeco(*cascadeVertices[0]) = std::move(trk_px);
2514 trk_pyDeco(*cascadeVertices[0]) = std::move(trk_py);
2515 trk_pzDeco(*cascadeVertices[0]) = std::move(trk_pz);
2516 trk_px_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Px();
2517 trk_py_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Py();
2518 trk_pz_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Pz();
2519
2520 result.push_back( std::make_pair(fit_result.release(),nullptr) );
2521 }
2522 }
2523 } // m_extraTrk1MassHypo<=0
2524 else { // m_extraTrk1MassHypo>0
2525 std::vector<double> massesExtra{m_extraTrk1MassHypo};
2526 std::vector<double> massesJXExtra = massesJX; massesJXExtra.push_back(m_extraTrk1MassHypo);
2527
2528 for(const xAOD::TrackParticle* tpExtra : *trackContainer) {
2529 if ( tpExtra->pt()<m_extraTrk1MinPt ) continue;
2530 if ( !m_trkSelector->decision(*tpExtra, nullptr) ) continue;
2531 // Check identical tracks in input
2532 if(std::find(tracksJX.cbegin(),tracksJX.cend(),tpExtra) != tracksJX.cend()) continue;
2533 if(std::find(tracksV0.cbegin(),tracksV0.cend(),tpExtra) != tracksV0.cend()) continue;
2534 if(tpExtra == disVtx.track) continue;
2535
2536 TLorentzVector tmp;
2537 tmp.SetPtEtaPhiM(tpExtra->pt(),tpExtra->eta(),tpExtra->phi(),m_extraTrk1MassHypo);
2538 double main_mass = (p4_moth+tmp).M();
2539 if(m_useImprovedMass) {
2540 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - (p4_moth - p4_disV).M() + m_massJpsi;
2541 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - (p4_moth - p4_disV).M() + m_massJX;
2542 if(m_massDisV>0) main_mass += - p4_disV.M() + m_massDisV;
2543 }
2544 if (main_mass < m_MassLower || main_mass > m_MassUpper) continue;
2545
2546 std::vector<const xAOD::TrackParticle*> tracksExtra{tpExtra};
2547 std::vector<const xAOD::TrackParticle*> tracksJXExtra = tracksJX; tracksJXExtra.push_back(tpExtra);
2548
2549 // Apply the user's settings to the fitter
2550 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
2551 // Robustness: http://cdsweb.cern.ch/record/685551
2552 int robustness = 0;
2553 m_iVertexFitter->setRobustness(robustness, *state);
2554 // Build up the topology
2555 // Vertex list
2556 std::vector<Trk::VertexID> vrtList;
2557 std::vector<Trk::VertexID> vrtList2;
2558 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
2559 // V0 vertex
2560 Trk::VertexID vID1;
2561 if (m_constrV0) {
2562 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,disVtx.V0type==KS?m_massKs:m_massLd);
2563 } else {
2564 vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
2565 }
2566 vrtList.push_back(vID1);
2567 // Displaced vertex
2568 Trk::VertexID vID2;
2569 if (m_constrDisV) {
2570 vID2 = m_iVertexFitter->nextVertex(tracks3,massesDis3,vrtList,*state,m_massDisV);
2571 } else {
2572 vID2 = m_iVertexFitter->nextVertex(tracks3,massesDis3,vrtList,*state);
2573 }
2574 vrtList2.push_back(vID2);
2575 Trk::VertexID vID3;
2576 if(m_JXSubVtx) {
2577 // JXExtra vertex
2578 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state);
2579 vrtList2.push_back(vID3);
2580 // Mother vertex
2581 if(m_constrMainV) {
2582 m_iVertexFitter->nextVertex(tracksExtra,massesExtra,vrtList2,*state,m_massMainV);
2583 } else {
2584 m_iVertexFitter->nextVertex(tracksExtra,massesExtra,vrtList2,*state);
2585 }
2586 }
2587 else { // m_JXSubVtx=false
2588 // Mother vertex includes just one subvertex (DisV) and JX tracks + extra track
2589 if(m_constrMainV) {
2590 vID3 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList2,*state,m_massMainV);
2591 } else {
2592 vID3 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList2,*state);
2593 }
2594 }
2595 if (m_constrJX && m_jxDaug_num>2) {
2596 std::vector<Trk::VertexID> cnstV;
2597 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
2598 ATH_MSG_WARNING("addMassConstraint for JX failed");
2599 }
2600 }
2601 if (m_constrJpsi) {
2602 std::vector<Trk::VertexID> cnstV;
2603 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
2604 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
2605 }
2606 }
2607 if (m_constrX && m_jxDaug_num==4) {
2608 std::vector<Trk::VertexID> cnstV;
2609 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksX,cnstV,*state,m_massX).isSuccess() ) {
2610 ATH_MSG_WARNING("addMassConstraint for X failed");
2611 }
2612 }
2613 // Do the work
2614 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) );
2615
2616 if (fit_result) {
2617 for(auto& v : fit_result->vertices()) {
2618 if(v->nTrackParticles()==0) {
2619 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
2620 v->setTrackParticleLinks(nullLinkVector);
2621 }
2622 }
2623 // reset links to original tracks
2624 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
2625
2626 // necessary to prevent memory leak
2627 fit_result->setSVOwnership(true);
2628
2629 // Chi2/DOF cut
2630 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
2631 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
2632
2633 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
2634 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
2635 size_t iMoth = cascadeVertices.size()-1;
2636 double lxy_SV1_sub = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
2637 double lxy_SV1 = m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[iMoth]);
2638
2639 if(chi2CutPassed && lxy_SV1>m_lxyDisV_cut && lxy_SV1_sub>m_lxyV0_cut) {
2640 chi2_V0_decor(*cascadeVertices[0]) = disVtx.V0vtx->chiSquared();
2641 ndof_V0_decor(*cascadeVertices[0]) = disVtx.V0vtx->numberDoF();
2642 if(disVtx.V0type==LAMBDA) {
2643 type_V0_decor(*cascadeVertices[0]) = "Lambda";
2644 }
2645 else if(disVtx.V0type==LAMBDABAR) {
2646 type_V0_decor(*cascadeVertices[0]) = "Lambdabar";
2647 }
2648 else if(disVtx.V0type==KS) {
2649 type_V0_decor(*cascadeVertices[0]) = "Ks";
2650 }
2651 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*disVtx.V0vtx) ? mAcc_gfit(*disVtx.V0vtx) : 0;
2652 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*disVtx.V0vtx) ? mAcc_gmass(*disVtx.V0vtx) : -1;
2653 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*disVtx.V0vtx) ? mAcc_gmasserr(*disVtx.V0vtx) : -1;
2654 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*disVtx.V0vtx) ? mAcc_gchisq(*disVtx.V0vtx) : 999999;
2655 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*disVtx.V0vtx) ? mAcc_gndof(*disVtx.V0vtx) : 0;
2656 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*disVtx.V0vtx) ? mAcc_gprob(*disVtx.V0vtx) : -1;
2657 trk_px.clear(); trk_py.clear(); trk_pz.clear();
2658 trk_px.push_back( disVtx.p4_V0track1.Px() ); trk_px.push_back( disVtx.p4_V0track2.Px() );
2659 trk_py.push_back( disVtx.p4_V0track1.Py() ); trk_py.push_back( disVtx.p4_V0track2.Py() );
2660 trk_pz.push_back( disVtx.p4_V0track1.Pz() ); trk_pz.push_back( disVtx.p4_V0track2.Pz() );
2661 trk_pxDeco(*cascadeVertices[0]) = trk_px;
2662 trk_pyDeco(*cascadeVertices[0]) = trk_py;
2663 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
2664 trk_px_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Px();
2665 trk_py_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Py();
2666 trk_pz_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Pz();
2667
2668 result.push_back( std::make_pair(fit_result.release(),nullptr) );
2669 }
2670 }
2671 } // loop over trackContainer
2672 } // m_extraTrk1MassHypo>0
2673
2674 if(pv_xAOD) {
2675 for(auto cascade_info_pair : result) {
2676 if(cascade_info_pair.first && cascade_info_pair.first->getParticleMoms().size()>0) {
2677 size_t index = cascade_info_pair.first->getParticleMoms().size() - 1;
2678 const std::vector<TLorentzVector> &mom = cascade_info_pair.first->getParticleMoms()[index];
2679 const Amg::MatrixX &cov = cascade_info_pair.first->getCovariance()[index];
2680 const xAOD::Vertex* mainVertex = cascade_info_pair.first->vertices()[index];
2681 xAOD::BPhysHypoHelper vtx(m_hypoName, mainVertex);
2682 bool isInDefaultPVCont = false;
2683 for(const xAOD::Vertex* pvVtx : *defaultPVContainer) {
2684 if(pv_xAOD == pvVtx) { isInDefaultPVCont = true; break; }
2685 }
2686 if(isInDefaultPVCont) vtx.setPv( pv_xAOD, defaultPVContainer, pvtype );
2687 else vtx.setPv( pv_xAOD, pvContainer, pvtype );
2688 if(origPv_xAOD) vtx.setOrigPv( origPv_xAOD, defaultPVContainer, pvtype );
2689 vtx.setLxy ( m_CascadeTools->lxy (mom, vtx.vtx(), pv_xAOD), pvtype );
2690 vtx.setLxyErr ( m_CascadeTools->lxyError (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2691 vtx.setA0 ( m_CascadeTools->a0 (mom, vtx.vtx(), pv_xAOD), pvtype );
2692 vtx.setA0Err ( m_CascadeTools->a0Error (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2693 vtx.setA0xy ( m_CascadeTools->a0xy (mom, vtx.vtx(), pv_xAOD), pvtype );
2694 vtx.setA0xyErr( m_CascadeTools->a0xyError(mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2695 vtx.setZ0 ( m_CascadeTools->a0z (mom, vtx.vtx(), pv_xAOD), pvtype );
2696 vtx.setZ0Err ( m_CascadeTools->a0zError (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
2697 vtx.setRefitPVStatus( 0, pvtype );
2698 // Proper decay times
2699 vtx.setTau( m_CascadeTools->tau(mom, vtx.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
2700 vtx.setTauErr( m_CascadeTools->tauError(mom, cov, vtx.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
2701 vtx.setTau( m_CascadeTools->tau(mom, vtx.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
2702 vtx.setTauErr( m_CascadeTools->tauError(mom, cov, vtx.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
2703 }
2704 }
2705 }
2706
2707 return result;
2708 }
2709
2710 void JpsiXPlusDisplaced::fitV0Container(const EventContext& ctx, xAOD::VertexContainer* V0ContainerNew, const std::vector<const xAOD::TrackParticle*>& selectedTracks, const std::vector<const xAOD::TrackParticleContainer*>& trackCols) const {
2711
2712 SG::Decorator<std::string> mDec_type("Type_V0Vtx");
2713 SG::Decorator<int> mDec_gfit("gamma_fit");
2714 SG::Decorator<float> mDec_gmass("gamma_mass");
2715 SG::Decorator<float> mDec_gmasserr("gamma_massError");
2716 SG::Decorator<float> mDec_gchisq("gamma_chisq");
2717 SG::Decorator<int> mDec_gndof("gamma_ndof");
2718 SG::Decorator<float> mDec_gprob("gamma_probability");
2719
2720 std::vector<const xAOD::TrackParticle*> posTracks;
2721 std::vector<const xAOD::TrackParticle*> negTracks;
2722 for(const xAOD::TrackParticle* TP : selectedTracks) {
2723 if(TP->charge()>0) posTracks.push_back(TP);
2724 else negTracks.push_back(TP);
2725 }
2726
2727 for(const xAOD::TrackParticle* TP1 : posTracks) {
2728 const Trk::Perigee& aPerigee1 = TP1->perigeeParameters();
2729 for(const xAOD::TrackParticle* TP2 : negTracks) {
2730 const Trk::Perigee& aPerigee2 = TP2->perigeeParameters();
2731 int sflag(0), errorcode(0);
2732 Amg::Vector3D startingPoint = m_vertexEstimator->getCirclesIntersectionPoint(&aPerigee1,&aPerigee2,sflag,errorcode);
2733 if (errorcode != 0) {startingPoint(0) = 0.0; startingPoint(1) = 0.0; startingPoint(2) = 0.0;}
2734
2735 if (errorcode == 0 || errorcode == 5 || errorcode == 6 || errorcode == 8) {
2736 Trk::PerigeeSurface perigeeSurface(startingPoint);
2737 const Trk::TrackParameters* extrapolatedPerigee1 = m_extrapolator->extrapolate(ctx,TP1->perigeeParameters(), perigeeSurface).release();
2738 const Trk::TrackParameters* extrapolatedPerigee2 = m_extrapolator->extrapolate(ctx,TP2->perigeeParameters(), perigeeSurface).release();
2739 std::vector<std::unique_ptr<const Trk::TrackParameters> > cleanup;
2740 if(!extrapolatedPerigee1) extrapolatedPerigee1 = &TP1->perigeeParameters();
2741 else cleanup.push_back(std::unique_ptr<const Trk::TrackParameters>(extrapolatedPerigee1));
2742 if(!extrapolatedPerigee2) extrapolatedPerigee2 = &TP2->perigeeParameters();
2743 else cleanup.push_back(std::unique_ptr<const Trk::TrackParameters>(extrapolatedPerigee2));
2744 if(extrapolatedPerigee1 && extrapolatedPerigee2) {
2745 bool pass = false;
2746 TLorentzVector v1; TLorentzVector v2;
2747 if(!pass) {
2748 v1.SetXYZM(extrapolatedPerigee1->momentum().x(),extrapolatedPerigee1->momentum().y(),extrapolatedPerigee1->momentum().z(),m_mass_proton);
2749 v2.SetXYZM(extrapolatedPerigee2->momentum().x(),extrapolatedPerigee2->momentum().y(),extrapolatedPerigee2->momentum().z(),m_mass_pion);
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_proton);
2755 if((v1+v2).M()>1030.0 && (v1+v2).M()<1200.0) pass = true;
2756 }
2757 if(!pass) {
2758 v1.SetXYZM(extrapolatedPerigee1->momentum().x(),extrapolatedPerigee1->momentum().y(),extrapolatedPerigee1->momentum().z(),m_mass_pion);
2759 v2.SetXYZM(extrapolatedPerigee2->momentum().x(),extrapolatedPerigee2->momentum().y(),extrapolatedPerigee2->momentum().z(),m_mass_pion);
2760 if((v1+v2).M()>430.0 && (v1+v2).M()<565.0) pass = true;
2761 }
2762 if(pass) {
2763 std::vector<const xAOD::TrackParticle*> tracksV0;
2764 tracksV0.push_back(TP1); tracksV0.push_back(TP2);
2765 std::unique_ptr<xAOD::Vertex> V0vtx = m_iV0Fitter->fit(ctx, tracksV0, startingPoint);
2766 if(V0vtx && V0vtx->chiSquared()>=0) {
2767 double chi2DOF = V0vtx->chiSquared()/V0vtx->numberDoF();
2768 if(chi2DOF>m_chi2cut_V0) continue;
2769
2770 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);
2771 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);
2772 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);
2773 if(massSig_V0_Lambda1<=massSig_V0_Lambda2 && massSig_V0_Lambda1<=massSig_V0_Ks) {
2774 mDec_type(*V0vtx.get()) = "Lambda";
2775 }
2776 else if(massSig_V0_Lambda2<=massSig_V0_Lambda1 && massSig_V0_Lambda2<=massSig_V0_Ks) {
2777 mDec_type(*V0vtx.get()) = "Lambdabar";
2778 }
2779 else if(massSig_V0_Ks<=massSig_V0_Lambda1 && massSig_V0_Ks<=massSig_V0_Lambda2) {
2780 mDec_type(*V0vtx.get()) = "Ks";
2781 }
2782
2783 int gamma_fit = 0; int gamma_ndof = 0; double gamma_chisq = 999999.;
2784 double gamma_prob = -1., gamma_mass = -1., gamma_massErr = -1.;
2785 std::unique_ptr<xAOD::Vertex> gammaVtx = m_iGammaFitter->fit(ctx, tracksV0, m_V0Tools->vtx(V0vtx.get()));
2786 if (gammaVtx) {
2787 gamma_fit = 1;
2788 gamma_mass = m_V0Tools->invariantMass(gammaVtx.get(),m_mass_e,m_mass_e);
2789 gamma_massErr = m_V0Tools->invariantMassError(gammaVtx.get(),m_mass_e,m_mass_e);
2790 gamma_chisq = m_V0Tools->chisq(gammaVtx.get());
2791 gamma_ndof = m_V0Tools->ndof(gammaVtx.get());
2792 gamma_prob = m_V0Tools->vertexProbability(gammaVtx.get());
2793 }
2794 mDec_gfit(*V0vtx.get()) = gamma_fit;
2795 mDec_gmass(*V0vtx.get()) = gamma_mass;
2796 mDec_gmasserr(*V0vtx.get()) = gamma_massErr;
2797 mDec_gchisq(*V0vtx.get()) = gamma_chisq;
2798 mDec_gndof(*V0vtx.get()) = gamma_ndof;
2799 mDec_gprob(*V0vtx.get()) = gamma_prob;
2800
2801 xAOD::BPhysHelper V0_helper(V0vtx.get());
2802 V0_helper.setRefTrks(); // AOD only method
2803
2804 if(not trackCols.empty()){
2805 try {
2806 JpsiUpsilonCommon::RelinkVertexTracks(trackCols, V0vtx.get());
2807 } catch (std::runtime_error const& e) {
2808 ATH_MSG_ERROR(e.what());
2809 return;
2810 }
2811 }
2812
2813 V0ContainerNew->push_back(std::move(V0vtx));
2814 }
2815 }
2816 }
2817 }
2818 }
2819 }
2820 }
2821
2822 template<size_t NTracks>
2823 const xAOD::Vertex* JpsiXPlusDisplaced::FindVertex(const xAOD::VertexContainer* cont, const xAOD::Vertex* v) const {
2824 for (const xAOD::Vertex* v1 : *cont) {
2825 assert(v1->nTrackParticles() == NTracks);
2826 std::array<const xAOD::TrackParticle*, NTracks> a1;
2827 std::array<const xAOD::TrackParticle*, NTracks> a2;
2828 for(size_t i=0; i<NTracks; i++){
2829 a1[i] = v1->trackParticle(i);
2830 a2[i] = v->trackParticle(i);
2831 }
2832 std::sort(a1.begin(), a1.end());
2833 std::sort(a2.begin(), a2.end());
2834 if(a1 == a2) return v1;
2835 }
2836 return nullptr;
2837 }
2838}
ATH_CHECK
#define ATH_CHECK
Evaluate an expression and check for errors.
Definition AthCheckMacros.h:40
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x,...)
Definition AthMsgStreamMacros.h:43
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x,...)
Definition AthMsgStreamMacros.h:47
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x,...)
Definition AthMsgStreamMacros.h:46
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x,...)
Definition AthMsgStreamMacros.h:48
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:1283
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:1314
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:1247
DerivationFramework::JpsiXPlusDisplaced::FindVertex
const xAOD::Vertex * FindVertex(const xAOD::VertexContainer *cont, const xAOD::Vertex *v) const
Definition JpsiXPlusDisplaced.cxx:2823
DerivationFramework::JpsiXPlusDisplaced::getXiCandidate
XiCandidate getXiCandidate(const EventContext &ctx, const xAOD::Vertex *V0vtx, const V0Enum V0, const xAOD::TrackParticle *track3) const
Definition JpsiXPlusDisplaced.cxx:1180
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:2710
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:1170
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:1230
DerivationFramework::JpsiXPlusDisplaced::m_V0Tools
PublicToolHandle< Trk::V0Tools > m_V0Tools
Definition JpsiXPlusDisplaced.h:185
SG::Accessor
Helper class to provide type-safe access to aux data.
Definition Control/AthContainers/AthContainers/Accessor.h:68
SG::ConstAccessor::isAvailable
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
SG::ReadHandleKey
Property holding a SG store/key/clid from which a ReadHandle is made.
Definition StoreGate/StoreGate/ReadHandleKey.h:40
SG::ReadHandle::ptr
const_pointer_type ptr()
Dereference the pointer.
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
SG::ReadHandle::cptr
const_pointer_type cptr()
Dereference the pointer.
SG::WriteHandleKey
Property holding a SG store/key/clid from which a WriteHandle is made.
Definition StoreGate/StoreGate/WriteHandleKey.h:40
SG::WriteHandle::cptr
const_pointer_type cptr() const
Dereference the pointer.
SG::WriteHandle::record
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
SG::WriteHandle::ptr
pointer_type ptr()
Dereference the pointer.
Trk::ParametersBase::momentum
const Amg::Vector3D & momentum() const
Access method for the momentum.
Trk::PerigeeSurface
Class describing the Line to which the Perigee refers to.
Definition PerigeeSurface.h:43
xAOD::BPhysHelper::setZ0
float setZ0(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
longitudinal impact parameter
Definition BPhysHelper.cxx:1067
xAOD::BPhysHelper::origPv
const xAOD::Vertex * origPv(const pv_type vertexType=BPhysHelper::PV_MIN_A0)
original PV
Definition BPhysHelper.cxx:807
xAOD::BPhysHelper::setLxyErr
bool setLxyErr(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
its error
Definition BPhysHelper.cxx:899
xAOD::BPhysHelper::setLxy
bool setLxy(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
Set the transverse decay distance and its error measured between the refitted primary vertex of type ...
Definition BPhysHelper.cxx:888
xAOD::BPhysHelper::setZ0Err
float setZ0Err(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
longitudinal impact parameter error
Definition BPhysHelper.cxx:1078
xAOD::BPhysHelper::vtx
const xAOD::Vertex * vtx() const
Getter method for the cached xAOD::Vertex.
Definition BPhysHelper.h:108
xAOD::BPhysHelper::setA0
float setA0(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
Set the 3D and transverse impact parameters and their error.
Definition BPhysHelper.cxx:1023
xAOD::BPhysHelper::refTrks
const std::vector< TVector3 > & refTrks()
Returns refitted track momenta.
Definition BPhysHelper.cxx:142
xAOD::BPhysHelper::setPv
bool setPv(const xAOD::Vertex *pv, const xAOD::VertexContainer *vertexContainer, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
Set the refitted collision vertex of type pv_type.
Definition BPhysHelper.cxx:831
xAOD::BPhysHelper::pv
const xAOD::Vertex * pv(const pv_type vertexType=BPhysHelper::PV_MIN_A0)
Get the refitted collision vertex of type pv_type.
Definition BPhysHelper.cxx:796
xAOD::BPhysHelper::refTrk
TVector3 refTrk(const size_t index)
Returns i-th refitted track 3-momentum.
Definition BPhysHelper.cxx:126
xAOD::BPhysHelper::pv_type
pv_type
: Enum type of the PV
Definition BPhysHelper.h:475
xAOD::BPhysHelper::setRefTrks
bool setRefTrks(std::vector< float > px, std::vector< float > py, std::vector< float > pz)
Sets refitted track momenta.
Definition BPhysHelper.cxx:286
xAOD::BPhysHelper::setA0Err
float setA0Err(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
3D impact parameter error
Definition BPhysHelper.cxx:1034
xAOD::BPhysHelper::setRefitPVStatus
bool setRefitPVStatus(int code, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
Set the exitCode of the refitter for vertex of type pv_type.
Definition BPhysHelper.cxx:844
xAOD::BPhysHelper::setA0xyErr
float setA0xyErr(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
transverse impact parameter error
Definition BPhysHelper.cxx:1056
xAOD::BPhysHelper::nRefTrks
int nRefTrks()
Returns number of stored refitted track momenta.
Definition BPhysHelper.cxx:115
xAOD::BPhysHelper::setA0xy
float setA0xy(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
transverse impact parameter
Definition BPhysHelper.cxx:1045
xAOD::BPhysHelper::setOrigPv
bool setOrigPv(const xAOD::Vertex *pv, const xAOD::VertexContainer *vertexContainer, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
Set the original collision vertex of type pv_type.
Definition BPhysHelper.cxx:818
xAOD::BPhysHypoHelper::setMass
bool setMass(const float val)
Set given invariant mass and its error.
Definition BPhysHypoHelper.cxx:52
xAOD::BPhysHypoHelper::setPass
bool setPass(bool passVal)
get the pass flag for this hypothesis
Definition BPhysHypoHelper.cxx:367
xAOD::BPhysHypoHelper::setTau
bool setTau(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0, const tau_type tauType=BPhysHypoHelper::TAU_CONST_MASS)
: Set the proper decay time and error.
Definition BPhysHypoHelper.cxx:143
xAOD::BPhysHypoHelper::setMassErr
bool setMassErr(const float val)
invariant mass error
Definition BPhysHypoHelper.cxx:57
xAOD::BPhysHypoHelper::setTauErr
bool setTauErr(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0, const tau_type tauType=BPhysHypoHelper::TAU_CONST_MASS)
proper decay time error
Definition BPhysHypoHelper.cxx:179
xAOD::TrackParticle_v1::perigeeParameters
const Trk::Perigee & perigeeParameters() const
Returns the Trk::MeasuredPerigee track parameters.
Definition TrackParticle_v1.cxx: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:48
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
SG::Decorator
SG::Decorator< T, ALLOC > Decorator
Helper class to provide type-safe access to aux data, specialized for JaggedVecElt.
Definition AuxElement.h:576
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:264
xAOD::numberOfPixelHits
@ numberOfPixelHits
these are the pixel hits, including the b-layer [unit8_t].
Definition TrackingPrimitives.h:255
type
Generated on for ATLAS Offline Software by doxygen 1.17.0