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