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