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(const std::string& 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) {
713  ATH_MSG_ERROR("CascadeInfo is null");
714  continue;
715  }
716 
717  const std::vector<xAOD::Vertex*> &cascadeVertices = cascade_info->vertices();
718  if(cascadeVertices.size() != topoN) ATH_MSG_ERROR("Incorrect number of vertices");
719  for(size_t i=0; i<topoN; i++) {
720  if(cascadeVertices[i]==nullptr) ATH_MSG_ERROR("Error null vertex");
721  }
722 
723  cascade_info->setSVOwnership(false); // Prevent Container from deleting vertices
724  const auto mainVertex = cascadeVertices[topoN-1]; // this is the mother vertex
725  const std::vector< std::vector<TLorentzVector> > &moms = cascade_info->getParticleMoms();
726 
727  // Identify the input JX
728  int ijx = m_JXSubVtx ? topoN-2 : topoN-1;
729  if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0 && m_extraTrk3MassHypo>0) ijx = topoN-1;
730  const xAOD::Vertex* jxVtx(nullptr);
731  if(m_jxDaug_num==4) jxVtx = FindVertex<4>(jxContainer.ptr(), cascadeVertices[ijx]);
732  else if(m_jxDaug_num==3) jxVtx = FindVertex<3>(jxContainer.ptr(), cascadeVertices[ijx]);
733  else jxVtx = FindVertex<2>(jxContainer.ptr(), cascadeVertices[ijx]);
734 
735  xAOD::BPhysHypoHelper vtx(m_hypoName, mainVertex);
736 
737  // Get refitted track momenta from all vertices, charged tracks only
738  BPhysPVCascadeTools::SetVectorInfo(vtx, cascade_info);
739  vtx.setPass(true);
740 
741  //
742  // Decorate main vertex
743  //
744  // mass, mass error
745  // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/VxCascadeInfo.h
746  BPHYS_CHECK( vtx.setMass(m_CascadeTools->invariantMass(moms[topoN-1])) );
747  BPHYS_CHECK( vtx.setMassErr(m_CascadeTools->invariantMassError(moms[topoN-1],cascade_info->getCovariance()[topoN-1])) );
748  // pt and pT error (the default pt of mainVertex is != the pt of the full cascade fit!)
749  Pt_decor(*mainVertex) = m_CascadeTools->pT(moms[topoN-1]);
750  PtErr_decor(*mainVertex) = m_CascadeTools->pTError(moms[topoN-1],cascade_info->getCovariance()[topoN-1]);
751  // chi2 and ndof (the default chi2 of mainVertex is != the chi2 of the full cascade fit!)
752  chi2_decor(*mainVertex) = cascade_info->fitChi2();
753  ndof_decor(*mainVertex) = cascade_info->nDoF();
754 
755  if(m_disVDaug_num==2) {
756  // decorate the newly fitted V0 vertex
757  lxy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],mainVertex);
758  lxyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxyError(moms[0],cascade_info->getCovariance()[0],cascadeVertices[0],mainVertex);
759  a0z_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0z(moms[0],cascadeVertices[0],mainVertex);
760  a0zErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0zError(moms[0],cascade_info->getCovariance()[0],cascadeVertices[0],mainVertex);
761  a0xy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xy(moms[0],cascadeVertices[0],mainVertex);
762  a0xyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xyError(moms[0],cascade_info->getCovariance()[0],cascadeVertices[0],mainVertex);
763  }
764  else {
765  // decorate the newly fitted V0 vertex
766  lxy_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
767  lxyErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->lxyError(moms[0],cascade_info->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
768  a0z_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0z(moms[0],cascadeVertices[0],cascadeVertices[1]);
769  a0zErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0zError(moms[0],cascade_info->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
770  a0xy_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0xy(moms[0],cascadeVertices[0],cascadeVertices[1]);
771  a0xyErr_SV0_decor(*cascadeVertices[0]) = m_CascadeTools->a0xyError(moms[0],cascade_info->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
772 
773  // decorate the newly fitted disV vertex
774  lxy_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->lxy(moms[1],cascadeVertices[1],mainVertex);
775  lxyErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->lxyError(moms[1],cascade_info->getCovariance()[1],cascadeVertices[1],mainVertex);
776  a0xy_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0z(moms[1],cascadeVertices[1],mainVertex);
777  a0xyErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0zError(moms[1],cascade_info->getCovariance()[1],cascadeVertices[1],mainVertex);
778  a0z_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0xy(moms[1],cascadeVertices[1],mainVertex);
779  a0zErr_SV1_decor(*cascadeVertices[1]) = m_CascadeTools->a0xyError(moms[1],cascade_info->getCovariance()[1],cascadeVertices[1],mainVertex);
780  }
781 
782  if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0 && m_extraTrk3MassHypo>0) {
783  lxy_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->lxy(moms[1],cascadeVertices[1],mainVertex);
784  lxyErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->lxyError(moms[1],cascade_info->getCovariance()[1],cascadeVertices[1],mainVertex);
785  a0z_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0z(moms[1],cascadeVertices[1],mainVertex);
786  a0zErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0zError(moms[1],cascade_info->getCovariance()[1],cascadeVertices[1],mainVertex);
787  a0xy_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0xy(moms[1],cascadeVertices[1],mainVertex);
788  a0xyErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0xyError(moms[1],cascade_info->getCovariance()[1],cascadeVertices[1],mainVertex);
789  }
790  // decorate the newly fitted JX vertex
791  else if(m_JXSubVtx) {
792  lxy_SV2_decor(*cascadeVertices[ijx]) = m_CascadeTools->lxy(moms[ijx],cascadeVertices[ijx],mainVertex);
793  lxyErr_SV2_decor(*cascadeVertices[ijx]) = m_CascadeTools->lxyError(moms[ijx],cascade_info->getCovariance()[ijx],cascadeVertices[ijx],mainVertex);
794  a0z_SV2_decor(*cascadeVertices[ijx]) = m_CascadeTools->a0z(moms[ijx],cascadeVertices[ijx],mainVertex);
795  a0zErr_SV2_decor(*cascadeVertices[ijx]) = m_CascadeTools->a0zError(moms[ijx],cascade_info->getCovariance()[ijx],cascadeVertices[ijx],mainVertex);
796  a0xy_SV2_decor(*cascadeVertices[ijx]) = m_CascadeTools->a0xy(moms[ijx],cascadeVertices[ijx],mainVertex);
797  a0xyErr_SV2_decor(*cascadeVertices[ijx]) = m_CascadeTools->a0xyError(moms[ijx],cascade_info->getCovariance()[ijx],cascadeVertices[ijx],mainVertex);
798  }
799 
800  chi2_V2_decor(*cascadeVertices[ijx]) = m_V0Tools->chisq(jxVtx);
801  ndof_V2_decor(*cascadeVertices[ijx]) = m_V0Tools->ndof(jxVtx);
802 
803  double Mass_Moth = m_CascadeTools->invariantMass(moms[topoN-1]);
804  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));
805 
806  for(size_t i=0; i<topoN; i++) {
807  VtxWriteHandles[i].ptr()->push_back(cascadeVertices[i]);
808  }
809 
810  // Set links to cascade vertices
811  VertexLinkVector precedingVertexLinks;
812  VertexLink vertexLink1;
813  vertexLink1.setElement(cascadeVertices[0]);
814  vertexLink1.setStorableObject(*VtxWriteHandles[0].ptr());
815  if( vertexLink1.isValid() ) precedingVertexLinks.push_back( vertexLink1 );
816  if(topoN>=3) {
817  VertexLink vertexLink2;
818  vertexLink2.setElement(cascadeVertices[1]);
819  vertexLink2.setStorableObject(*VtxWriteHandles[1].ptr());
820  if( vertexLink2.isValid() ) precedingVertexLinks.push_back( vertexLink2 );
821  }
822  if(topoN==4) {
823  VertexLink vertexLink3;
824  vertexLink3.setElement(cascadeVertices[2]);
825  vertexLink3.setStorableObject(*VtxWriteHandles[2].ptr());
826  if( vertexLink3.isValid() ) precedingVertexLinks.push_back( vertexLink3 );
827  }
828  CascadeLinksDecor(*mainVertex) = precedingVertexLinks;
829  } // loop over cascadeinfoContainer
830 
831  // Deleting cascadeinfo since this won't be stored.
832  for (auto cascade_info : cascadeinfoContainer) delete cascade_info;
833 
834  return StatusCode::SUCCESS;
835  }
836 
837  bool JpsiXPlusDisplaced::d0Pass(const xAOD::TrackParticle* track, const xAOD::Vertex* PV) const {
838  bool pass = false;
839  const EventContext& ctx = Gaudi::Hive::currentContext();
840  std::unique_ptr<Trk::Perigee> per = m_trackToVertexTool->perigeeAtVertex(ctx, *track, PV->position());
841  if(!per) return pass;
842  double d0 = per->parameters()[Trk::d0];
843  double sig_d0 = sqrt((*per->covariance())(0,0));
844  if(std::abs(d0/sig_d0) > m_d0_cut) pass = true;
845  return pass;
846  }
847 
848  JpsiXPlusDisplaced::XiCandidate JpsiXPlusDisplaced::getXiCandidate(const xAOD::Vertex* V0vtx, const V0Enum V0, const xAOD::TrackParticle* track3) const {
849  XiCandidate disVtx;
850  // Check overlap
851  std::vector<const xAOD::TrackParticle*> tracksV0;
852  tracksV0.reserve(V0vtx->nTrackParticles());
853  for(size_t i=0; i<V0vtx->nTrackParticles(); i++) tracksV0.push_back(V0vtx->trackParticle(i));
854  if(std::find(tracksV0.cbegin(), tracksV0.cend(), track3) != tracksV0.cend()) return disVtx;
855 
856  std::vector<double> massesV0;
857  if(V0==LAMBDA) {
858  massesV0 = m_massesV0_ppi;
859  }
860  else if(V0==LAMBDABAR) {
861  massesV0 = m_massesV0_pip;
862  }
863  else if(V0==KS) {
864  massesV0 = m_massesV0_pipi;
865  }
866  xAOD::BPhysHelper V0_helper(V0vtx);
867  TLorentzVector p4_v0, tmp;
868  for(int i=0; i<V0_helper.nRefTrks(); i++) p4_v0 += V0_helper.refTrk(i,massesV0[i]);
869  tmp.SetPtEtaPhiM(track3->pt(),track3->eta(),track3->phi(),m_disVDaug3MassHypo);
870  // A rough mass window cut, as V0 and the track are not from a common vertex
871  if((p4_v0+tmp).M() < m_DisplacedMassLower-120. || (p4_v0+tmp).M() > m_DisplacedMassUpper+120.) return disVtx;
872 
873  std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
874  int robustness = 0;
875  m_iVertexFitter->setRobustness(robustness, *state);
876  std::vector<Trk::VertexID> vrtList;
877  // V0 vertex
878  Trk::VertexID vID = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
879  vrtList.push_back(vID);
880  // Mother vertex
881  std::vector<const xAOD::TrackParticle*> tracksDis3{track3};
882  std::vector<double> massesDis3{m_disVDaug3MassHypo};
883  m_iVertexFitter->nextVertex(tracksDis3,massesDis3,vrtList,*state);
884  // Do the work
885  std::unique_ptr<Trk::VxCascadeInfo> cascade_info = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state) );
886  if(cascade_info) {
887  cascade_info->setSVOwnership(true);
888  double chi2NDF = cascade_info->fitChi2()/cascade_info->nDoF();
889  if(m_chi2cut_DisV<=0 || chi2NDF < m_chi2cut_DisV) {
890  const std::vector<std::vector<TLorentzVector> > &moms = cascade_info->getParticleMoms();
891  const std::vector<xAOD::Vertex*> &cascadeVertices = cascade_info->vertices();
892  double lxy_SV1_sub = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
893  if(lxy_SV1_sub > m_lxyV0_cut) {
894  TLorentzVector totalMom;
895  for(size_t it=0; it<moms[1].size(); it++) totalMom += moms[1][it];
896  if(totalMom.M()>m_DisplacedMassLower && totalMom.M()<m_DisplacedMassUpper) {
897  disVtx.track = track3; disVtx.V0vtx = V0vtx; disVtx.V0type = V0;
898  disVtx.chi2NDF = chi2NDF; disVtx.p4_V0track1 = moms[0][0];
899  disVtx.p4_V0track2 = moms[0][1]; disVtx.p4_disVtrack = moms[1][0];
900  }
901  }
902  }
903  }
904  return disVtx;
905  }
906 
907  JpsiXPlusDisplaced::MesonCandidate JpsiXPlusDisplaced::getEtacCandidate(const xAOD::Vertex* V0vtx, const V0Enum V0, const xAOD::TrackParticle* extraTrk1, const xAOD::TrackParticle* extraTrk2) const {
908  MesonCandidate etac;
909  // Check overlap
910  std::vector<const xAOD::TrackParticle*> tracksV0;
911  tracksV0.reserve(V0vtx->nTrackParticles());
912  for(size_t i=0; i<V0vtx->nTrackParticles(); i++) tracksV0.push_back(V0vtx->trackParticle(i));
913  if(std::find(tracksV0.cbegin(), tracksV0.cend(), extraTrk1) != tracksV0.cend()) return etac;
914  if(std::find(tracksV0.cbegin(), tracksV0.cend(), extraTrk2) != tracksV0.cend()) return etac;
915 
916  std::vector<double> massesV0;
917  if(V0==LAMBDA) {
918  massesV0 = m_massesV0_ppi;
919  }
920  else if(V0==LAMBDABAR) {
921  massesV0 = m_massesV0_pip;
922  }
923  else if(V0==KS) {
924  massesV0 = m_massesV0_pipi;
925  }
926  xAOD::BPhysHelper V0_helper(V0vtx);
927  TLorentzVector p4_v0, tmp1, tmp2;
928  for(int i=0; i<V0_helper.nRefTrks(); i++) p4_v0 += V0_helper.refTrk(i,massesV0[i]);
929  tmp1.SetPtEtaPhiM(extraTrk1->pt(),extraTrk1->eta(),extraTrk1->phi(),m_extraTrk1MassHypo);
930  tmp2.SetPtEtaPhiM(extraTrk2->pt(),extraTrk2->eta(),extraTrk2->phi(),m_extraTrk2MassHypo);
931  if((p4_v0+tmp1+tmp2).M() < m_V0ExtraMassLower || (p4_v0+tmp1+tmp2).M() > m_V0ExtraMassUpper) return etac;
932 
933  std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
934  int robustness = 0;
935  m_iVertexFitter->setRobustness(robustness, *state);
936  std::vector<Trk::VertexID> vrtList;
937  // V0 vertex
938  Trk::VertexID vID = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
939  vrtList.push_back(vID);
940  // Mother vertex
941  std::vector<const xAOD::TrackParticle*> extraTracks{extraTrk1, extraTrk2};
942  std::vector<double> extraMasses{m_extraTrk1MassHypo, m_extraTrk2MassHypo};
943  m_iVertexFitter->nextVertex(extraTracks,extraMasses,vrtList,*state);
944  // Do the work
945  std::unique_ptr<Trk::VxCascadeInfo> cascade_info = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state) );
946  if(cascade_info) {
947  cascade_info->setSVOwnership(true);
948  double chi2NDF = cascade_info->fitChi2()/cascade_info->nDoF();
949  if(m_chi2cut_V0Extra<=0 || chi2NDF < m_chi2cut_V0Extra) {
950  const std::vector<std::vector<TLorentzVector> > &moms = cascade_info->getParticleMoms();
951  const std::vector<xAOD::Vertex*> &cascadeVertices = cascade_info->vertices();
952  double lxy = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
953  if(lxy > m_lxyV0_cut) {
954  TLorentzVector totalMom;
955  for(size_t it=0; it<moms[1].size(); it++) totalMom += moms[1][it];
956  etac.V0vtx = V0vtx; etac.V0type = V0;
957  etac.extraTrack1 = extraTrk1; etac.extraTrack2 = extraTrk2;
958  etac.chi2NDF = chi2NDF; etac.pt = totalMom.Pt();
959  }
960  }
961  }
962  return etac;
963  }
964 
965  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 {
966  MesonCandidate Dpm;
967  // Check overlap
968  std::vector<const xAOD::TrackParticle*> tracksJX;
969  tracksJX.reserve(JXvtx->nTrackParticles());
970  for(size_t i=0; i<JXvtx->nTrackParticles(); i++) tracksJX.push_back(JXvtx->trackParticle(i));
971  if(std::find(tracksJX.cbegin(), tracksJX.cend(), extraTrk1) != tracksJX.cend()) return Dpm;
972  if(std::find(tracksJX.cbegin(), tracksJX.cend(), extraTrk2) != tracksJX.cend()) return Dpm;
973  if(std::find(tracksJX.cbegin(), tracksJX.cend(), extraTrk3) != tracksJX.cend()) return Dpm;
974 
975  TLorentzVector tmp1, tmp2, tmp3;
976  tmp1.SetPtEtaPhiM(extraTrk1->pt(),extraTrk1->eta(),extraTrk1->phi(),m_extraTrk1MassHypo);
977  tmp2.SetPtEtaPhiM(extraTrk2->pt(),extraTrk2->eta(),extraTrk2->phi(),m_extraTrk2MassHypo);
978  tmp3.SetPtEtaPhiM(extraTrk3->pt(),extraTrk3->eta(),extraTrk3->phi(),m_extraTrk3MassHypo);
979  if((tmp1+tmp2+tmp3).M() < m_DpmMassLower || (tmp1+tmp2+tmp3).M() > m_DpmMassUpper) return Dpm;
980 
981  std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
982  int robustness = 0;
983  m_iVertexFitter->setRobustness(robustness, *state);
984  std::vector<Trk::VertexID> vrtList;
985  // Dpm vertex
986  std::vector<const xAOD::TrackParticle*> extraTracks{extraTrk1, extraTrk2, extraTrk3};
987  std::vector<double> extraMasses{m_extraTrk1MassHypo, m_extraTrk2MassHypo, m_extraTrk3MassHypo};
988  Trk::VertexID vID = m_iVertexFitter->startVertex(extraTracks,extraMasses,*state);
989  vrtList.push_back(vID);
990  // Mother vertex
991  m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
992  // Do the work
993  std::unique_ptr<Trk::VxCascadeInfo> cascade_info = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state) );
994  if(cascade_info) {
995  cascade_info->setSVOwnership(true);
996  double chi2NDF = cascade_info->fitChi2()/cascade_info->nDoF();
997  if(m_chi2cut_JXDpm<=0 || chi2NDF < m_chi2cut_JXDpm) {
998  const std::vector<std::vector<TLorentzVector> > &moms = cascade_info->getParticleMoms();
999  const std::vector<xAOD::Vertex*> &cascadeVertices = cascade_info->vertices();
1000  double lxy = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
1001  if(lxy > m_lxyDpm_cut) {
1002  Dpm.extraTrack1 = extraTrk1; Dpm.extraTrack2 = extraTrk2; Dpm.extraTrack3 = extraTrk3;
1003  Dpm.chi2NDF = chi2NDF; Dpm.pt = moms[1][tracksJX.size()].Pt();
1004  }
1005  }
1006  }
1007  return Dpm;
1008  }
1009 
1010  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 {
1011  MesonCandidate Dstpm;
1012  // Check overlap
1013  std::vector<const xAOD::TrackParticle*> tracksJX;
1014  tracksJX.reserve(JXvtx->nTrackParticles());
1015  for(size_t i=0; i<JXvtx->nTrackParticles(); i++) tracksJX.push_back(JXvtx->trackParticle(i));
1016  if(std::find(tracksJX.cbegin(), tracksJX.cend(), extraTrk1) != tracksJX.cend()) return Dstpm;
1017  if(std::find(tracksJX.cbegin(), tracksJX.cend(), extraTrk2) != tracksJX.cend()) return Dstpm;
1018  if(std::find(tracksJX.cbegin(), tracksJX.cend(), extraTrk3) != tracksJX.cend()) return Dstpm;
1019 
1020  TLorentzVector tmp1, tmp2, tmp3;
1021  tmp1.SetPtEtaPhiM(extraTrk1->pt(),extraTrk1->eta(),extraTrk1->phi(),m_extraTrk1MassHypo);
1022  tmp2.SetPtEtaPhiM(extraTrk2->pt(),extraTrk2->eta(),extraTrk2->phi(),m_extraTrk2MassHypo);
1023  tmp3.SetPtEtaPhiM(extraTrk3->pt(),extraTrk3->eta(),extraTrk3->phi(),m_extraTrk3MassHypo);
1024  if((tmp1+tmp2).M() < m_D0MassLower || (tmp1+tmp2).M() > m_D0MassUpper) return Dstpm;
1025  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;
1026 
1027  std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
1028  int robustness = 0;
1029  m_iVertexFitter->setRobustness(robustness, *state);
1030  std::vector<Trk::VertexID> vrtList;
1031  // D0 vertex
1032  std::vector<const xAOD::TrackParticle*> extraTracks{extraTrk1, extraTrk2};
1033  std::vector<double> extraMasses{m_extraTrk1MassHypo, m_extraTrk2MassHypo};
1034  Trk::VertexID vID = m_iVertexFitter->startVertex(extraTracks,extraMasses,*state);
1035  vrtList.push_back(vID);
1036  // Mother vertex
1037  std::vector<const xAOD::TrackParticle*> tracksJXExtra = tracksJX;
1038  tracksJXExtra.push_back(extraTrk3);
1039  std::vector<double> massesJXExtra = massesJX;
1040  massesJXExtra.push_back(m_extraTrk3MassHypo);
1041  m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList,*state);
1042  // Do the work
1043  std::unique_ptr<Trk::VxCascadeInfo> cascade_info = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state) );
1044  if(cascade_info) {
1045  cascade_info->setSVOwnership(true);
1046  double chi2NDF = cascade_info->fitChi2()/cascade_info->nDoF();
1047  if(m_chi2cut_JXDstpm<=0 || chi2NDF < m_chi2cut_JXDstpm) {
1048  const std::vector<std::vector<TLorentzVector> > &moms = cascade_info->getParticleMoms();
1049  const std::vector<xAOD::Vertex*> &cascadeVertices = cascade_info->vertices();
1050  double lxy = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
1051  if(lxy > m_lxyD0_cut) {
1052  TLorentzVector totalMom;
1053  for(size_t it=tracksJX.size(); it<moms[1].size(); it++) totalMom += moms[1][it];
1054  Dstpm.extraTrack1 = extraTrk1; Dstpm.extraTrack2 = extraTrk2; Dstpm.extraTrack3 = extraTrk3;
1055  Dstpm.chi2NDF = chi2NDF; Dstpm.pt = totalMom.Pt();
1056  }
1057  }
1058  }
1059  return Dstpm;
1060  }
1061 
1062  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 {
1063  std::vector<Trk::VxCascadeInfo*> result;
1064 
1065  std::vector<const xAOD::TrackParticle*> tracksJX;
1066  tracksJX.reserve(JXvtx->nTrackParticles());
1067  for(size_t i=0; i<JXvtx->nTrackParticles(); i++) tracksJX.push_back(JXvtx->trackParticle(i));
1068  if (tracksJX.size() != massesJX.size()) {
1069  ATH_MSG_ERROR("Problems with JX input: number of tracks or track mass inputs is not correct!");
1070  return result;
1071  }
1072  // Check identical tracks in input
1073  if(std::find(tracksJX.cbegin(), tracksJX.cend(), V0vtx->trackParticle(0)) != tracksJX.cend()) return result;
1074  if(std::find(tracksJX.cbegin(), tracksJX.cend(), V0vtx->trackParticle(1)) != tracksJX.cend()) return result;
1075  std::vector<const xAOD::TrackParticle*> tracksV0;
1076  tracksV0.reserve(V0vtx->nTrackParticles());
1077  for(size_t j=0; j<V0vtx->nTrackParticles(); j++) tracksV0.push_back(V0vtx->trackParticle(j));
1078 
1079  std::vector<const xAOD::TrackParticle*> tracksJpsi{tracksJX[0], tracksJX[1]};
1080  std::vector<const xAOD::TrackParticle*> tracksX;
1081  if(m_jxDaug_num>=3) tracksX.push_back(tracksJX[2]);
1082  if(m_jxDaug_num==4) tracksX.push_back(tracksJX[3]);
1083 
1084  std::vector<double> massesV0;
1085  if(V0==LAMBDA) {
1086  massesV0 = m_massesV0_ppi;
1087  }
1088  else if(V0==LAMBDABAR) {
1089  massesV0 = m_massesV0_pip;
1090  }
1091  else if(V0==KS) {
1092  massesV0 = m_massesV0_pipi;
1093  }
1094 
1095  TLorentzVector p4_moth, p4_v0, tmp;
1096  for(size_t it=0; it<JXvtx->nTrackParticles(); it++) {
1097  tmp.SetPtEtaPhiM(JXvtx->trackParticle(it)->pt(), JXvtx->trackParticle(it)->eta(), JXvtx->trackParticle(it)->phi(), massesJX[it]);
1098  p4_moth += tmp;
1099  }
1100  xAOD::BPhysHelper V0_helper(V0vtx);
1101  for(int it=0; it<V0_helper.nRefTrks(); it++) {
1102  p4_moth += V0_helper.refTrk(it,massesV0[it]);
1103  p4_v0 += V0_helper.refTrk(it,massesV0[it]);
1104  }
1105 
1106  SG::AuxElement::Decorator<float> chi2_V1_decor("ChiSquared_V1");
1107  SG::AuxElement::Decorator<int> ndof_V1_decor("nDoF_V1");
1108  SG::AuxElement::Decorator<std::string> type_V1_decor("Type_V1");
1109 
1110  SG::AuxElement::Accessor<int> mAcc_gfit("gamma_fit");
1111  SG::AuxElement::Accessor<float> mAcc_gmass("gamma_mass");
1112  SG::AuxElement::Accessor<float> mAcc_gmasserr("gamma_massError");
1113  SG::AuxElement::Accessor<float> mAcc_gchisq("gamma_chisq");
1114  SG::AuxElement::Accessor<int> mAcc_gndof("gamma_ndof");
1115  SG::AuxElement::Accessor<float> mAcc_gprob("gamma_probability");
1116 
1117  SG::AuxElement::Decorator<int> mDec_gfit("gamma_fit");
1118  SG::AuxElement::Decorator<float> mDec_gmass("gamma_mass");
1119  SG::AuxElement::Decorator<float> mDec_gmasserr("gamma_massError");
1120  SG::AuxElement::Decorator<float> mDec_gchisq("gamma_chisq");
1121  SG::AuxElement::Decorator<int> mDec_gndof("gamma_ndof");
1122  SG::AuxElement::Decorator<float> mDec_gprob("gamma_probability");
1123  SG::AuxElement::Decorator< std::vector<float> > trk_pxDeco("TrackPx_V0nc");
1124  SG::AuxElement::Decorator< std::vector<float> > trk_pyDeco("TrackPy_V0nc");
1125  SG::AuxElement::Decorator< std::vector<float> > trk_pzDeco("TrackPz_V0nc");
1126 
1127  std::vector<float> trk_px;
1128  std::vector<float> trk_py;
1129  std::vector<float> trk_pz;
1130 
1131  if(m_extraTrk1MassHypo<=0) {
1132  if (p4_moth.M() < m_MassLower || p4_moth.M() > m_MassUpper) return result;
1133 
1134  // Apply the user's settings to the fitter
1135  std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
1136  // Robustness: http://cdsweb.cern.ch/record/685551
1137  int robustness = 0;
1138  m_iVertexFitter->setRobustness(robustness, *state);
1139  // Build up the topology
1140  // Vertex list
1141  std::vector<Trk::VertexID> vrtList;
1142  // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
1143  // V0 vertex
1144  Trk::VertexID vID1;
1145  if (m_constrV0) {
1146  vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,m_massV0);
1147  } else {
1148  vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1149  }
1150  vrtList.push_back(vID1);
1151  Trk::VertexID vID2;
1152  if(m_JXSubVtx) {
1153  // JX vertex
1154  if (m_constrJX && m_jxDaug_num>2) {
1155  vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state,m_massJX);
1156  } else {
1157  vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state);
1158  }
1159  vrtList.push_back(vID2);
1160  // Mother vertex including JX and V0
1161  std::vector<const xAOD::TrackParticle*> tp;
1162  std::vector<double> tp_masses;
1163  if(m_constrMainV) {
1164  m_iVertexFitter->nextVertex(tp,tp_masses,vrtList,*state,m_massMainV);
1165  } else {
1166  m_iVertexFitter->nextVertex(tp,tp_masses,vrtList,*state);
1167  }
1168  }
1169  else { // m_JXSubVtx=false
1170  // Mother vertex including JX and V0
1171  if(m_constrMainV) {
1172  vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state,m_massMainV);
1173  } else {
1174  vID2 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
1175  }
1176  if (m_constrJX && m_jxDaug_num>2) {
1177  std::vector<Trk::VertexID> cnstV;
1178  if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
1179  ATH_MSG_WARNING("addMassConstraint for JX failed");
1180  }
1181  }
1182  }
1183  if (m_constrJpsi) {
1184  std::vector<Trk::VertexID> cnstV;
1185  if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
1186  ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
1187  }
1188  }
1189  if (m_constrX && m_jxDaug_num==4 && m_massX>0) {
1190  std::vector<Trk::VertexID> cnstV;
1191  if ( !m_iVertexFitter->addMassConstraint(vID2,tracksX,cnstV,*state,m_massX).isSuccess() ) {
1192  ATH_MSG_WARNING("addMassConstraint for X failed");
1193  }
1194  }
1195  // Do the work
1196  std::unique_ptr<Trk::VxCascadeInfo> fit_result = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state) );
1197 
1198  if (fit_result) {
1199  for(auto& v : fit_result->vertices()) {
1200  if(v->nTrackParticles()==0) {
1201  std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
1202  v->setTrackParticleLinks(nullLinkVector);
1203  }
1204  }
1205  // reset links to original tracks
1206  BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
1207 
1208  // necessary to prevent memory leak
1209  fit_result->setSVOwnership(true);
1210 
1211  // Chi2/DOF cut
1212  double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
1213  bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
1214 
1215  const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
1216  const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
1217  size_t iMoth = cascadeVertices.size()-1;
1218  double lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[iMoth]);
1219  if(chi2CutPassed && lxy_SV1>m_lxyV0_cut) {
1220  chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
1221  ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
1222  if(V0==LAMBDA) {
1223  type_V1_decor(*cascadeVertices[0]) = "Lambda";
1224  }
1225  else if(V0==LAMBDABAR) {
1226  type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
1227  }
1228  else if(V0==KS) {
1229  type_V1_decor(*cascadeVertices[0]) = "Ks";
1230  }
1231  mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
1232  mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
1233  mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
1234  mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
1235  mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
1236  mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
1237  trk_px.clear(); trk_py.clear(); trk_pz.clear();
1238  trk_px.reserve(V0_helper.nRefTrks());
1239  trk_py.reserve(V0_helper.nRefTrks());
1240  trk_pz.reserve(V0_helper.nRefTrks());
1241  for(auto&& vec3 : V0_helper.refTrks()) {
1242  trk_px.push_back( vec3.Px() );
1243  trk_py.push_back( vec3.Py() );
1244  trk_pz.push_back( vec3.Pz() );
1245  }
1246  trk_pxDeco(*cascadeVertices[0]) = trk_px;
1247  trk_pyDeco(*cascadeVertices[0]) = trk_py;
1248  trk_pzDeco(*cascadeVertices[0]) = trk_pz;
1249 
1250  result.push_back( fit_result.release() );
1251  }
1252  }
1253  } // m_extraTrk1MassHypo<=0
1254  else if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo<=0) {
1255  std::vector<double> massesJXExtra = massesJX;
1256  massesJXExtra.push_back(m_extraTrk1MassHypo);
1257 
1258  for(const xAOD::TrackParticle* tpExtra : *trackContainer) {
1259  if( tpExtra->pt()<m_extraTrk1MinPt ) continue;
1260  if( !m_trkSelector->decision(*tpExtra, nullptr) ) continue;
1261  // Check identical tracks in input
1262  if(std::find(tracksJX.cbegin(),tracksJX.cend(),tpExtra) != tracksJX.cend()) continue;
1263  if(std::find(tracksV0.cbegin(),tracksV0.cend(),tpExtra) != tracksV0.cend()) continue;
1264 
1265  TLorentzVector tmp;
1266  tmp.SetPtEtaPhiM(tpExtra->pt(),tpExtra->eta(),tpExtra->phi(),m_extraTrk1MassHypo);
1267  if((p4_moth+tmp).M() < m_MassLower || (p4_moth+tmp).M() > m_MassUpper) continue;
1268 
1269  std::vector<const xAOD::TrackParticle*> tracksJXExtra = tracksJX;
1270  tracksJXExtra.push_back(tpExtra);
1271 
1272  // Apply the user's settings to the fitter
1273  std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
1274  // Robustness: http://cdsweb.cern.ch/record/685551
1275  int robustness = 0;
1276  m_iVertexFitter->setRobustness(robustness, *state);
1277  // Build up the topology
1278  // Vertex list
1279  std::vector<Trk::VertexID> vrtList;
1280  // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
1281  // V0 vertex
1282  Trk::VertexID vID1;
1283  if (m_constrV0) {
1284  vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,m_massV0);
1285  } else {
1286  vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1287  }
1288  vrtList.push_back(vID1);
1289  Trk::VertexID vID2;
1290  if(m_JXSubVtx) {
1291  // JXExtra vertex
1292  vID2 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,*state);
1293  vrtList.push_back(vID2);
1294  // Mother vertex includes two subvertices: V0, JX+extra track
1295  std::vector<const xAOD::TrackParticle*> tp;
1296  std::vector<double> tp_masses;
1297  if(m_constrMainV) {
1298  m_iVertexFitter->nextVertex(tp,tp_masses,vrtList,*state,m_massMainV);
1299  } else {
1300  m_iVertexFitter->nextVertex(tp,tp_masses,vrtList,*state);
1301  }
1302  }
1303  else { // m_JXSubVtx=false
1304  // Mother vertex includes one subvertex (V0) and JX tracks + extra track
1305  if(m_constrMainV) {
1306  vID2 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList,*state,m_massMainV);
1307  } else {
1308  vID2 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList,*state);
1309  }
1310  }
1311  if (m_constrJX && m_jxDaug_num>2) {
1312  std::vector<Trk::VertexID> cnstV;
1313  if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
1314  ATH_MSG_WARNING("addMassConstraint for JX failed");
1315  }
1316  }
1317  if (m_constrJpsi) {
1318  std::vector<Trk::VertexID> cnstV;
1319  if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
1320  ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
1321  }
1322  }
1323  if (m_constrX && m_jxDaug_num==4 && m_massX>0) {
1324  std::vector<Trk::VertexID> cnstV;
1325  if ( !m_iVertexFitter->addMassConstraint(vID2,tracksX,cnstV,*state,m_massX).isSuccess() ) {
1326  ATH_MSG_WARNING("addMassConstraint for X failed");
1327  }
1328  }
1329  // Do the work
1330  std::unique_ptr<Trk::VxCascadeInfo> fit_result = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state) );
1331 
1332  if (fit_result) {
1333  for(auto& v : fit_result->vertices()) {
1334  if(v->nTrackParticles()==0) {
1335  std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
1336  v->setTrackParticleLinks(nullLinkVector);
1337  }
1338  }
1339  // reset links to original tracks
1340  BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
1341 
1342  // necessary to prevent memory leak
1343  fit_result->setSVOwnership(true);
1344 
1345  // Chi2/DOF cut
1346  double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
1347  bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
1348  const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
1349  const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
1350  size_t iMoth = cascadeVertices.size()-1;
1351  double lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[iMoth]);
1352  if(chi2CutPassed && lxy_SV1>m_lxyV0_cut) {
1353  chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
1354  ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
1355  if(V0==LAMBDA) {
1356  type_V1_decor(*cascadeVertices[0]) = "Lambda";
1357  }
1358  else if(V0==LAMBDABAR) {
1359  type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
1360  }
1361  else if(V0==KS) {
1362  type_V1_decor(*cascadeVertices[0]) = "Ks";
1363  }
1364  mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
1365  mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
1366  mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
1367  mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
1368  mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
1369  mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
1370  trk_px.clear(); trk_py.clear(); trk_pz.clear();
1371  trk_px.reserve(V0_helper.nRefTrks());
1372  trk_py.reserve(V0_helper.nRefTrks());
1373  trk_pz.reserve(V0_helper.nRefTrks());
1374  for(auto&& vec3 : V0_helper.refTrks()) {
1375  trk_px.push_back( vec3.Px() );
1376  trk_py.push_back( vec3.Py() );
1377  trk_pz.push_back( vec3.Pz() );
1378  }
1379  trk_pxDeco(*cascadeVertices[0]) = trk_px;
1380  trk_pyDeco(*cascadeVertices[0]) = trk_py;
1381  trk_pzDeco(*cascadeVertices[0]) = trk_pz;
1382 
1383  result.push_back( fit_result.release() );
1384  }
1385  }
1386  } // loop over trackContainer
1387  } // m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo<=0
1388  else if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0 && m_extraTrk3MassHypo<=0) {
1389  std::vector<const xAOD::TrackParticle*> tracksPlus;
1390  std::vector<const xAOD::TrackParticle*> tracksMinus;
1391  for(const xAOD::TrackParticle* tpExtra : *trackContainer) {
1392  if( tpExtra->pt() < std::fmin(m_extraTrk1MinPt,m_extraTrk2MinPt) ) continue;
1393  if( !m_trkSelector->decision(*tpExtra, nullptr) ) continue;
1394  // Check identical tracks in input
1395  if(std::find(tracksJX.cbegin(),tracksJX.cend(),tpExtra) != tracksJX.cend()) continue;
1396  if(std::find(tracksV0.cbegin(),tracksV0.cend(),tpExtra) != tracksV0.cend()) continue;
1397  if(tpExtra->charge()>0) {
1398  tracksPlus.push_back(tpExtra);
1399  }
1400  else {
1401  tracksMinus.push_back(tpExtra);
1402  }
1403  }
1404 
1405  MesonCandidateVector etacCandidates(m_maxMesonCandidates, m_MesonPtOrdering);
1406  TLorentzVector p4_ExtraTrk1, p4_ExtraTrk2;
1407  for(const xAOD::TrackParticle* tp1 : tracksPlus) {
1408  for(const xAOD::TrackParticle* tp2 : tracksMinus) {
1409  if((tp1->pt()>m_extraTrk1MinPt && tp2->pt()>m_extraTrk2MinPt) ||
1410  (tp1->pt()>m_extraTrk2MinPt && tp2->pt()>m_extraTrk1MinPt)) {
1411  p4_ExtraTrk1.SetPtEtaPhiM(tp1->pt(), tp1->eta(), tp1->phi(), m_extraTrk1MassHypo);
1412  p4_ExtraTrk2.SetPtEtaPhiM(tp2->pt(), tp2->eta(), tp2->phi(), m_extraTrk2MassHypo);
1413  if((p4_moth+p4_ExtraTrk1+p4_ExtraTrk2).M() < m_MassLower || (p4_moth+p4_ExtraTrk1+p4_ExtraTrk2).M() > m_MassUpper) continue;
1414  auto etac = getEtacCandidate(V0vtx,V0,tp1,tp2);
1415  if(etac.V0vtx) etacCandidates.push_back(etac);
1416  }
1417  }
1418  }
1419 
1420  std::vector<double> massesJXExtra = massesJX;
1421  massesJXExtra.push_back(m_extraTrk1MassHypo);
1422  massesJXExtra.push_back(m_extraTrk2MassHypo);
1423 
1424  for(auto&& etac : etacCandidates.vector()) {
1425  std::vector<const xAOD::TrackParticle*> tracksExtra{etac.extraTrack1, etac.extraTrack2};
1426  std::vector<const xAOD::TrackParticle*> tracksJXExtra = tracksJX;
1427  tracksJXExtra.push_back(etac.extraTrack1);
1428  tracksJXExtra.push_back(etac.extraTrack2);
1429 
1430  // Apply the user's settings to the fitter
1431  std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
1432  // Robustness: http://cdsweb.cern.ch/record/685551
1433  int robustness = 0;
1434  m_iVertexFitter->setRobustness(robustness, *state);
1435  // Build up the topology
1436  // Vertex list
1437  std::vector<Trk::VertexID> vrtList;
1438  // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
1439  // V0 vertex
1440  Trk::VertexID vID1;
1441  if (m_constrV0) {
1442  vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,m_massV0);
1443  } else {
1444  vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1445  }
1446  vrtList.push_back(vID1);
1447  Trk::VertexID vID2;
1448  if(m_JXSubVtx) {
1449  // JXExtra vertex
1450  vID2 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,*state);
1451  vrtList.push_back(vID2);
1452  // Mother vertex includes two subvertices: V0, JX+extra tracks
1453  std::vector<const xAOD::TrackParticle*> tp;
1454  std::vector<double> tp_masses;
1455  if(m_constrMainV) {
1456  m_iVertexFitter->nextVertex(tp,tp_masses,vrtList,*state,m_massMainV);
1457  } else {
1458  m_iVertexFitter->nextVertex(tp,tp_masses,vrtList,*state);
1459  }
1460  }
1461  else { // m_JXSubVtx=false
1462  // Mother vertex includes one subvertex (V0) and JX tracks + extra track
1463  if(m_constrMainV) {
1464  vID2 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList,*state,m_massMainV);
1465  } else {
1466  vID2 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList,*state);
1467  }
1468  }
1469  if (m_constrJX && m_jxDaug_num>2) {
1470  std::vector<Trk::VertexID> cnstV;
1471  if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
1472  ATH_MSG_WARNING("addMassConstraint for JX failed");
1473  }
1474  }
1475  if (m_constrJpsi) {
1476  std::vector<Trk::VertexID> cnstV;
1477  if ( !m_iVertexFitter->addMassConstraint(vID2,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
1478  ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
1479  }
1480  }
1481  if (m_constrX && m_jxDaug_num==4 && m_massX>0) {
1482  std::vector<Trk::VertexID> cnstV;
1483  if ( !m_iVertexFitter->addMassConstraint(vID2,tracksX,cnstV,*state,m_massX).isSuccess() ) {
1484  ATH_MSG_WARNING("addMassConstraint for X failed");
1485  }
1486  }
1487  if (m_constrV0Extra && m_massV0Extra>0) {
1488  std::vector<Trk::VertexID> cnstV{vID1};
1489  if ( !m_iVertexFitter->addMassConstraint(vID2,tracksExtra,cnstV,*state,m_massV0Extra).isSuccess() ) {
1490  ATH_MSG_WARNING("addMassConstraint for V0+extraTracks failed");
1491  }
1492  }
1493  // Do the work
1494  std::unique_ptr<Trk::VxCascadeInfo> fit_result = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state) );
1495 
1496  if (fit_result) {
1497  for(auto& v : fit_result->vertices()) {
1498  if(v->nTrackParticles()==0) {
1499  std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
1500  v->setTrackParticleLinks(nullLinkVector);
1501  }
1502  }
1503  // reset links to original tracks
1504  BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
1505 
1506  // necessary to prevent memory leak
1507  fit_result->setSVOwnership(true);
1508 
1509  // Chi2/DOF cut
1510  double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
1511  bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
1512  const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
1513  const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
1514  size_t iMoth = cascadeVertices.size()-1;
1515  double lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[iMoth]);
1516  if(chi2CutPassed && lxy_SV1>m_lxyV0_cut) {
1517  chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
1518  ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
1519  if(V0==LAMBDA) {
1520  type_V1_decor(*cascadeVertices[0]) = "Lambda";
1521  }
1522  else if(V0==LAMBDABAR) {
1523  type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
1524  }
1525  else if(V0==KS) {
1526  type_V1_decor(*cascadeVertices[0]) = "Ks";
1527  }
1528  mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
1529  mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
1530  mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
1531  mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
1532  mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
1533  mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
1534  trk_px.clear(); trk_py.clear(); trk_pz.clear();
1535  trk_px.reserve(V0_helper.nRefTrks());
1536  trk_py.reserve(V0_helper.nRefTrks());
1537  trk_pz.reserve(V0_helper.nRefTrks());
1538  for(auto&& vec3 : V0_helper.refTrks()) {
1539  trk_px.push_back( vec3.Px() );
1540  trk_py.push_back( vec3.Py() );
1541  trk_pz.push_back( vec3.Pz() );
1542  }
1543  trk_pxDeco(*cascadeVertices[0]) = trk_px;
1544  trk_pyDeco(*cascadeVertices[0]) = trk_py;
1545  trk_pzDeco(*cascadeVertices[0]) = trk_pz;
1546 
1547  result.push_back( fit_result.release() );
1548  }
1549  }
1550  } // loop over etacCandidates
1551  } // m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0 && m_extraTrk3MassHypo<=0
1552  else if(m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0 && m_extraTrk3MassHypo>0) {
1553  std::vector<const xAOD::TrackParticle*> tracksPlus;
1554  std::vector<const xAOD::TrackParticle*> tracksMinus;
1555  double minTrkPt = std::fmin(std::fmin(m_extraTrk1MinPt,m_extraTrk2MinPt),m_extraTrk3MinPt);
1556  for(const xAOD::TrackParticle* tpExtra : *trackContainer) {
1557  if( tpExtra->pt() < minTrkPt ) continue;
1558  if( !m_trkSelector->decision(*tpExtra, nullptr) ) continue;
1559  // Check identical tracks in input
1560  if(std::find(tracksJX.cbegin(),tracksJX.cend(),tpExtra) != tracksJX.cend()) continue;
1561  if(std::find(tracksV0.cbegin(),tracksV0.cend(),tpExtra) != tracksV0.cend()) continue;
1562  if(tpExtra->charge()>0) {
1563  tracksPlus.push_back(tpExtra);
1564  }
1565  else {
1566  tracksMinus.push_back(tpExtra);
1567  }
1568  }
1569 
1570  MesonCandidateVector DpmCandidates(m_maxMesonCandidates, m_MesonPtOrdering);
1571  MesonCandidateVector DstpmCandidates(m_maxMesonCandidates, m_MesonPtOrdering);
1572  TLorentzVector p4_ExtraTrk1, p4_ExtraTrk2, p4_ExtraTrk3;
1573  // +,-,- combination (D*-/D- -> K+ pi- pi-)
1574  for(const xAOD::TrackParticle* tp1 : tracksPlus) {
1575  if( tp1->pt() < m_extraTrk1MinPt ) continue;
1576  for(auto tp2Itr=tracksMinus.cbegin(); tp2Itr!=tracksMinus.cend(); ++tp2Itr) {
1577  const xAOD::TrackParticle* tp2 = *tp2Itr;
1578  for(auto tp3Itr=tp2Itr+1; tp3Itr!=tracksMinus.cend(); ++tp3Itr) {
1579  const xAOD::TrackParticle* tp3 = *tp3Itr;
1580  if((tp2->pt()>m_extraTrk2MinPt && tp3->pt()>m_extraTrk3MinPt) ||
1581  (tp2->pt()>m_extraTrk3MinPt && tp3->pt()>m_extraTrk2MinPt)) {
1582  p4_ExtraTrk1.SetPtEtaPhiM(tp1->pt(), tp1->eta(), tp1->phi(), m_extraTrk1MassHypo);
1583  p4_ExtraTrk2.SetPtEtaPhiM(tp2->pt(), tp2->eta(), tp2->phi(), m_extraTrk2MassHypo);
1584  p4_ExtraTrk3.SetPtEtaPhiM(tp3->pt(), tp3->eta(), tp3->phi(), m_extraTrk3MassHypo);
1585  if((p4_moth+p4_ExtraTrk1+p4_ExtraTrk2+p4_ExtraTrk3).M() < m_MassLower || (p4_moth+p4_ExtraTrk1+p4_ExtraTrk2+p4_ExtraTrk3).M() > m_MassUpper) continue;
1586  auto Dpm = getDpmCandidate(JXvtx,massesJX,tp1,tp2,tp3);
1587  if(Dpm.extraTrack1) DpmCandidates.push_back(Dpm);
1588  auto Dstpm = getDstpmCandidate(JXvtx,massesJX,tp1,tp2,tp3);
1589  if(Dstpm.extraTrack1) DstpmCandidates.push_back(Dstpm);
1590  }
1591  }
1592  }
1593  }
1594  // -,+,+ combination (D*+/D+ -> K- pi+ pi+)
1595  for(const xAOD::TrackParticle* tp1 : tracksMinus) {
1596  if( tp1->pt() < m_extraTrk1MinPt ) continue;
1597  for(auto tp2Itr=tracksPlus.cbegin(); tp2Itr!=tracksPlus.cend(); ++tp2Itr) {
1598  const xAOD::TrackParticle* tp2 = *tp2Itr;
1599  for(auto tp3Itr=tp2Itr+1; tp3Itr!=tracksPlus.cend(); ++tp3Itr) {
1600  const xAOD::TrackParticle* tp3 = *tp3Itr;
1601  if((tp2->pt()>m_extraTrk2MinPt && tp3->pt()>m_extraTrk3MinPt) ||
1602  (tp2->pt()>m_extraTrk3MinPt && tp3->pt()>m_extraTrk2MinPt)) {
1603  p4_ExtraTrk1.SetPtEtaPhiM(tp1->pt(), tp1->eta(), tp1->phi(), m_extraTrk1MassHypo);
1604  p4_ExtraTrk2.SetPtEtaPhiM(tp2->pt(), tp2->eta(), tp2->phi(), m_extraTrk2MassHypo);
1605  p4_ExtraTrk3.SetPtEtaPhiM(tp3->pt(), tp3->eta(), tp3->phi(), m_extraTrk3MassHypo);
1606  if((p4_moth+p4_ExtraTrk1+p4_ExtraTrk2+p4_ExtraTrk3).M() < m_MassLower || (p4_moth+p4_ExtraTrk1+p4_ExtraTrk2+p4_ExtraTrk3).M() > m_MassUpper) continue;
1607  auto Dpm = getDpmCandidate(JXvtx,massesJX,tp1,tp2,tp3);
1608  if(Dpm.extraTrack1) DpmCandidates.push_back(Dpm);
1609  auto Dstpm = getDstpmCandidate(JXvtx,massesJX,tp1,tp2,tp3);
1610  if(Dstpm.extraTrack1) DstpmCandidates.push_back(Dstpm);
1611  }
1612  }
1613  }
1614  }
1615 
1616  std::vector<double> massesExtra{m_extraTrk1MassHypo,m_extraTrk2MassHypo,m_extraTrk3MassHypo};
1617 
1618  for(auto&& Dpm : DpmCandidates.vector()) {
1619  std::vector<const xAOD::TrackParticle*> tracksExtra{Dpm.extraTrack1,Dpm.extraTrack2,Dpm.extraTrack3};
1620 
1621  // Apply the user's settings to the fitter
1622  std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
1623  // Robustness: http://cdsweb.cern.ch/record/685551
1624  int robustness = 0;
1625  m_iVertexFitter->setRobustness(robustness, *state);
1626  // Build up the topology
1627  // Vertex list
1628  std::vector<Trk::VertexID> vrtList;
1629  // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
1630  // V0 vertex
1631  Trk::VertexID vID1;
1632  if (m_constrV0) {
1633  vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,m_massV0);
1634  } else {
1635  vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1636  }
1637  vrtList.push_back(vID1);
1638  // Dpm vertex
1639  Trk::VertexID vID2;
1640  if (m_constrDpm && m_mass_Dpm>0) {
1641  vID2 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state,m_mass_Dpm);
1642  } else {
1643  vID2 = m_iVertexFitter->nextVertex(tracksExtra,massesExtra,*state);
1644  }
1645  vrtList.push_back(vID2);
1646  // Mother vertex includes two subvertices: V0, Dpm and JX tracks
1647  Trk::VertexID vID3;
1648  if(m_constrMainV) {
1649  vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state,m_massMainV);
1650  } else {
1651  vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
1652  }
1653  if (m_constrJX && m_jxDaug_num>2) {
1654  std::vector<Trk::VertexID> cnstV;
1655  if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
1656  ATH_MSG_WARNING("addMassConstraint for JX failed");
1657  }
1658  }
1659  if (m_constrJpsi) {
1660  std::vector<Trk::VertexID> cnstV;
1661  if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
1662  ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
1663  }
1664  }
1665  if (m_constrX && m_jxDaug_num==4 && m_massX>0) {
1666  std::vector<Trk::VertexID> cnstV;
1667  if ( !m_iVertexFitter->addMassConstraint(vID3,tracksX,cnstV,*state,m_massX).isSuccess() ) {
1668  ATH_MSG_WARNING("addMassConstraint for X failed");
1669  }
1670  }
1671  // Do the work
1672  std::unique_ptr<Trk::VxCascadeInfo> fit_result = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state) );
1673 
1674  if (fit_result) {
1675  for(auto& v : fit_result->vertices()) {
1676  if(v->nTrackParticles()==0) {
1677  std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
1678  v->setTrackParticleLinks(nullLinkVector);
1679  }
1680  }
1681  // reset links to original tracks
1682  BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
1683 
1684  // necessary to prevent memory leak
1685  fit_result->setSVOwnership(true);
1686 
1687  // Chi2/DOF cut
1688  double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
1689  bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
1690  const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
1691  const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
1692  size_t iMoth = cascadeVertices.size()-1;
1693  double lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[iMoth]);
1694  double lxy_SV2 = m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[iMoth]);
1695  if(chi2CutPassed && lxy_SV1>m_lxyV0_cut && lxy_SV2>m_lxyDpm_cut) {
1696  chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
1697  ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
1698  if(V0==LAMBDA) {
1699  type_V1_decor(*cascadeVertices[0]) = "Lambda";
1700  }
1701  else if(V0==LAMBDABAR) {
1702  type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
1703  }
1704  else if(V0==KS) {
1705  type_V1_decor(*cascadeVertices[0]) = "Ks";
1706  }
1707  mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
1708  mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
1709  mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
1710  mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
1711  mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
1712  mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
1713  trk_px.clear(); trk_py.clear(); trk_pz.clear();
1714  trk_px.reserve(V0_helper.nRefTrks());
1715  trk_py.reserve(V0_helper.nRefTrks());
1716  trk_pz.reserve(V0_helper.nRefTrks());
1717  for(auto&& vec3 : V0_helper.refTrks()) {
1718  trk_px.push_back( vec3.Px() );
1719  trk_py.push_back( vec3.Py() );
1720  trk_pz.push_back( vec3.Pz() );
1721  }
1722  trk_pxDeco(*cascadeVertices[0]) = trk_px;
1723  trk_pyDeco(*cascadeVertices[0]) = trk_py;
1724  trk_pzDeco(*cascadeVertices[0]) = trk_pz;
1725 
1726  result.push_back( fit_result.release() );
1727  }
1728  }
1729  } // loop over DpmCandidates
1730 
1731  std::vector<double> massesJXExtra = massesJX;
1732  massesJXExtra.push_back(m_extraTrk3MassHypo);
1733  std::vector<double> massesExtra12{m_extraTrk1MassHypo,m_extraTrk2MassHypo};
1734 
1735  for(auto&& Dstpm : DstpmCandidates.vector()) {
1736  std::vector<const xAOD::TrackParticle*> tracksExtra12{Dstpm.extraTrack1,Dstpm.extraTrack2};
1737  std::vector<const xAOD::TrackParticle*> tracksJXExtra = tracksJX;
1738  tracksJXExtra.push_back(Dstpm.extraTrack3);
1739 
1740  // Apply the user's settings to the fitter
1741  std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
1742  // Robustness: http://cdsweb.cern.ch/record/685551
1743  int robustness = 0;
1744  m_iVertexFitter->setRobustness(robustness, *state);
1745  // Build up the topology
1746  // Vertex list
1747  std::vector<Trk::VertexID> vrtList;
1748  // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
1749  // V0 vertex
1750  Trk::VertexID vID1;
1751  if (m_constrV0) {
1752  vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,m_massV0);
1753  } else {
1754  vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1755  }
1756  vrtList.push_back(vID1);
1757  // D0 vertex
1758  Trk::VertexID vID2;
1759  if (m_constrD0 && m_mass_D0>0) {
1760  vID2 = m_iVertexFitter->nextVertex(tracksExtra12,massesExtra12,*state,m_mass_D0);
1761  } else {
1762  vID2 = m_iVertexFitter->nextVertex(tracksExtra12,massesExtra12,*state);
1763  }
1764  vrtList.push_back(vID2);
1765  // Mother vertex includes two subvertices: V0, D0 and JX+extraTrk3 tracks
1766  Trk::VertexID vID3;
1767  if(m_constrMainV) {
1768  vID3 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList,*state,m_massMainV);
1769  } else {
1770  vID3 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList,*state);
1771  }
1772  if (m_constrJX && m_jxDaug_num>2) {
1773  std::vector<Trk::VertexID> cnstV;
1774  if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
1775  ATH_MSG_WARNING("addMassConstraint for JX failed");
1776  }
1777  }
1778  if (m_constrJpsi) {
1779  std::vector<Trk::VertexID> cnstV;
1780  if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
1781  ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
1782  }
1783  }
1784  if (m_constrX && m_jxDaug_num==4 && m_massX>0) {
1785  std::vector<Trk::VertexID> cnstV;
1786  if ( !m_iVertexFitter->addMassConstraint(vID3,tracksX,cnstV,*state,m_massX).isSuccess() ) {
1787  ATH_MSG_WARNING("addMassConstraint for X failed");
1788  }
1789  }
1790  if (m_constrDstpm && m_mass_Dstpm>0) {
1791  std::vector<Trk::VertexID> cnstV{vID2};
1792  std::vector<const xAOD::TrackParticle*> tracksExtra3{Dstpm.extraTrack3};
1793  if ( !m_iVertexFitter->addMassConstraint(vID3,tracksExtra3,cnstV,*state,m_mass_Dstpm).isSuccess() ) {
1794  ATH_MSG_WARNING("addMassConstraint for V0+extraTracks failed");
1795  }
1796  }
1797  // Do the work
1798  std::unique_ptr<Trk::VxCascadeInfo> fit_result = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state) );
1799 
1800  if (fit_result) {
1801  for(auto& v : fit_result->vertices()) {
1802  if(v->nTrackParticles()==0) {
1803  std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
1804  v->setTrackParticleLinks(nullLinkVector);
1805  }
1806  }
1807  // reset links to original tracks
1808  BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
1809 
1810  // necessary to prevent memory leak
1811  fit_result->setSVOwnership(true);
1812 
1813  // Chi2/DOF cut
1814  double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
1815  bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
1816  const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
1817  const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
1818  size_t iMoth = cascadeVertices.size()-1;
1819  double lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[iMoth]);
1820  double lxy_SV2 = m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[iMoth]);
1821  if(chi2CutPassed && lxy_SV1>m_lxyV0_cut && lxy_SV2>m_lxyD0_cut) {
1822  chi2_V1_decor(*cascadeVertices[0]) = V0vtx->chiSquared();
1823  ndof_V1_decor(*cascadeVertices[0]) = V0vtx->numberDoF();
1824  if(V0==LAMBDA) {
1825  type_V1_decor(*cascadeVertices[0]) = "Lambda";
1826  }
1827  else if(V0==LAMBDABAR) {
1828  type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
1829  }
1830  else if(V0==KS) {
1831  type_V1_decor(*cascadeVertices[0]) = "Ks";
1832  }
1833  mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V0vtx) ? mAcc_gfit(*V0vtx) : 0;
1834  mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V0vtx) ? mAcc_gmass(*V0vtx) : -1;
1835  mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V0vtx) ? mAcc_gmasserr(*V0vtx) : -1;
1836  mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V0vtx) ? mAcc_gchisq(*V0vtx) : 999999;
1837  mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V0vtx) ? mAcc_gndof(*V0vtx) : 0;
1838  mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V0vtx) ? mAcc_gprob(*V0vtx) : -1;
1839  trk_px.clear(); trk_py.clear(); trk_pz.clear();
1840  trk_px.reserve(V0_helper.nRefTrks());
1841  trk_py.reserve(V0_helper.nRefTrks());
1842  trk_pz.reserve(V0_helper.nRefTrks());
1843  for(auto&& vec3 : V0_helper.refTrks()) {
1844  trk_px.push_back( vec3.Px() );
1845  trk_py.push_back( vec3.Py() );
1846  trk_pz.push_back( vec3.Pz() );
1847  }
1848  trk_pxDeco(*cascadeVertices[0]) = trk_px;
1849  trk_pyDeco(*cascadeVertices[0]) = trk_py;
1850  trk_pzDeco(*cascadeVertices[0]) = trk_pz;
1851 
1852  result.push_back( fit_result.release() );
1853  }
1854  }
1855  } // loop over DstpmCandidates
1856  } // m_extraTrk1MassHypo>0 && m_extraTrk2MassHypo>0 && m_extraTrk3MassHypo>0
1857 
1858  return result;
1859  }
1860 
1861  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 {
1862  std::vector<Trk::VxCascadeInfo*> result;
1863 
1864  std::vector<const xAOD::TrackParticle*> tracksJX;
1865  tracksJX.reserve(JXvtx->nTrackParticles());
1866  for(size_t i=0; i<JXvtx->nTrackParticles(); i++) tracksJX.push_back(JXvtx->trackParticle(i));
1867  if (tracksJX.size() != massesJX.size()) {
1868  ATH_MSG_ERROR("Problems with JX input: number of tracks or track mass inputs is not correct!");
1869  return result;
1870  }
1871  // Check identical tracks in input
1872  if(std::find(tracksJX.cbegin(), tracksJX.cend(), disVtx.V0vtx->trackParticle(0)) != tracksJX.cend()) return result;
1873  if(std::find(tracksJX.cbegin(), tracksJX.cend(), disVtx.V0vtx->trackParticle(1)) != tracksJX.cend()) return result;
1874  std::vector<const xAOD::TrackParticle*> tracksV0;
1875  tracksV0.reserve(disVtx.V0vtx->nTrackParticles());
1876  for(size_t j=0; j<disVtx.V0vtx->nTrackParticles(); j++) tracksV0.push_back(disVtx.V0vtx->trackParticle(j));
1877 
1878  if(std::find(tracksJX.cbegin(), tracksJX.cend(), disVtx.track) != tracksJX.cend()) return result;
1879  std::vector<const xAOD::TrackParticle*> tracks3{disVtx.track};
1880  std::vector<double> massesDis3{m_disVDaug3MassHypo};
1881 
1882  std::vector<const xAOD::TrackParticle*> tracksJpsi{tracksJX[0], tracksJX[1]};
1883  std::vector<const xAOD::TrackParticle*> tracksX;
1884  if(m_jxDaug_num>=3) tracksX.push_back(tracksJX[2]);
1885  if(m_jxDaug_num==4) tracksX.push_back(tracksJX[3]);
1886 
1887  std::vector<double> massesV0;
1888  if(disVtx.V0type==LAMBDA) {
1889  massesV0 = m_massesV0_ppi;
1890  }
1891  else if(disVtx.V0type==LAMBDABAR) {
1892  massesV0 = m_massesV0_pip;
1893  }
1894  else if(disVtx.V0type==KS) {
1895  massesV0 = m_massesV0_pipi;
1896  }
1897 
1898  std::vector<double> massesDisV = massesV0;
1899  massesDisV.push_back(m_disVDaug3MassHypo);
1900 
1901  TLorentzVector p4_moth, tmp;
1902  for(size_t it=0; it<JXvtx->nTrackParticles(); it++) {
1903  tmp.SetPtEtaPhiM(JXvtx->trackParticle(it)->pt(), JXvtx->trackParticle(it)->eta(), JXvtx->trackParticle(it)->phi(), massesJX[it]);
1904  p4_moth += tmp;
1905  }
1906  p4_moth += disVtx.p4_V0track1; p4_moth += disVtx.p4_V0track2; p4_moth += disVtx.p4_disVtrack;
1907 
1908  SG::AuxElement::Decorator<float> chi2_V0_decor("ChiSquared_V0");
1909  SG::AuxElement::Decorator<int> ndof_V0_decor("nDoF_V0");
1910  SG::AuxElement::Decorator<std::string> type_V0_decor("Type_V0");
1911 
1912  SG::AuxElement::Accessor<int> mAcc_gfit("gamma_fit");
1913  SG::AuxElement::Accessor<float> mAcc_gmass("gamma_mass");
1914  SG::AuxElement::Accessor<float> mAcc_gmasserr("gamma_massError");
1915  SG::AuxElement::Accessor<float> mAcc_gchisq("gamma_chisq");
1916  SG::AuxElement::Accessor<int> mAcc_gndof("gamma_ndof");
1917  SG::AuxElement::Accessor<float> mAcc_gprob("gamma_probability");
1918 
1919  SG::AuxElement::Decorator<int> mDec_gfit("gamma_fit");
1920  SG::AuxElement::Decorator<float> mDec_gmass("gamma_mass");
1921  SG::AuxElement::Decorator<float> mDec_gmasserr("gamma_massError");
1922  SG::AuxElement::Decorator<float> mDec_gchisq("gamma_chisq");
1923  SG::AuxElement::Decorator<int> mDec_gndof("gamma_ndof");
1924  SG::AuxElement::Decorator<float> mDec_gprob("gamma_probability");
1925  SG::AuxElement::Decorator< std::vector<float> > trk_pxDeco("TrackPx_V0nc");
1926  SG::AuxElement::Decorator< std::vector<float> > trk_pyDeco("TrackPy_V0nc");
1927  SG::AuxElement::Decorator< std::vector<float> > trk_pzDeco("TrackPz_V0nc");
1928  SG::AuxElement::Decorator<float> trk_px_deco("TrackPx_DisVnc");
1929  SG::AuxElement::Decorator<float> trk_py_deco("TrackPy_DisVnc");
1930  SG::AuxElement::Decorator<float> trk_pz_deco("TrackPz_DisVnc");
1931 
1932  std::vector<float> trk_px;
1933  std::vector<float> trk_py;
1934  std::vector<float> trk_pz;
1935 
1936  if(m_extraTrk1MassHypo<=0) {
1937  if (p4_moth.M() < m_MassLower || p4_moth.M() > m_MassUpper) return result;
1938 
1939  // Apply the user's settings to the fitter
1940  std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
1941  // Robustness: http://cdsweb.cern.ch/record/685551
1942  int robustness = 0;
1943  m_iVertexFitter->setRobustness(robustness, *state);
1944  // Build up the topology
1945  // Vertex list
1946  std::vector<Trk::VertexID> vrtList;
1947  std::vector<Trk::VertexID> vrtList2;
1948  // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
1949  // V0 vertex
1950  Trk::VertexID vID1;
1951  if (m_constrV0) {
1952  vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,m_massV0);
1953  } else {
1954  vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
1955  }
1956  vrtList.push_back(vID1);
1957  // Displaced vertex
1958  Trk::VertexID vID2;
1959  if (m_constrDisV) {
1960  vID2 = m_iVertexFitter->nextVertex(tracks3,massesDis3,vrtList,*state,m_massDisV);
1961  } else {
1962  vID2 = m_iVertexFitter->nextVertex(tracks3,massesDis3,vrtList,*state);
1963  }
1964  vrtList2.push_back(vID2);
1965  Trk::VertexID vID3;
1966  if(m_JXSubVtx) {
1967  // JX vertex
1968  if (m_constrJX && m_jxDaug_num>2) {
1969  vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state,m_massJX);
1970  } else {
1971  vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state);
1972  }
1973  vrtList2.push_back(vID3);
1974  // Mother vertex includes two subvertices: DisV and JX
1975  std::vector<const xAOD::TrackParticle*> tp;
1976  std::vector<double> tp_masses;
1977  if(m_constrMainV) {
1978  m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2,*state,m_massMainV);
1979  } else {
1980  m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2,*state);
1981  }
1982  }
1983  else { // m_JXSubVtx=false
1984  // Mother vertex includes just one subvertex (DisV) and JX tracks
1985  if(m_constrMainV) {
1986  vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList2,*state,m_massMainV);
1987  } else {
1988  vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList2,*state);
1989  }
1990  if (m_constrJX && m_jxDaug_num>2) {
1991  std::vector<Trk::VertexID> cnstV;
1992  if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
1993  ATH_MSG_WARNING("addMassConstraint for JX failed");
1994  }
1995  }
1996  }
1997  if (m_constrJpsi) {
1998  std::vector<Trk::VertexID> cnstV;
1999  if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
2000  ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
2001  }
2002  }
2003  if (m_constrX && m_jxDaug_num==4 && m_massX>0) {
2004  std::vector<Trk::VertexID> cnstV;
2005  if ( !m_iVertexFitter->addMassConstraint(vID3,tracksX,cnstV,*state,m_massX).isSuccess() ) {
2006  ATH_MSG_WARNING("addMassConstraint for X failed");
2007  }
2008  }
2009  // Do the work
2010  std::unique_ptr<Trk::VxCascadeInfo> fit_result = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state) );
2011 
2012  if (fit_result) {
2013  for(auto& v : fit_result->vertices()) {
2014  if(v->nTrackParticles()==0) {
2015  std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
2016  v->setTrackParticleLinks(nullLinkVector);
2017  }
2018  }
2019  // reset links to original tracks
2020  BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
2021 
2022  // necessary to prevent memory leak
2023  fit_result->setSVOwnership(true);
2024 
2025  // Chi2/DOF cut
2026  double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
2027  bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
2028 
2029  const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
2030  const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
2031  size_t iMoth = cascadeVertices.size()-1;
2032  double lxy_SV1_sub = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
2033  double lxy_SV1 = m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[iMoth]);
2034 
2035  if(chi2CutPassed && lxy_SV1>m_lxyDisV_cut && lxy_SV1_sub>m_lxyV0_cut) {
2036  chi2_V0_decor(*cascadeVertices[0]) = disVtx.V0vtx->chiSquared();
2037  ndof_V0_decor(*cascadeVertices[0]) = disVtx.V0vtx->numberDoF();
2038  if(disVtx.V0type==LAMBDA) {
2039  type_V0_decor(*cascadeVertices[0]) = "Lambda";
2040  }
2041  else if(disVtx.V0type==LAMBDABAR) {
2042  type_V0_decor(*cascadeVertices[0]) = "Lambdabar";
2043  }
2044  else if(disVtx.V0type==KS) {
2045  type_V0_decor(*cascadeVertices[0]) = "Ks";
2046  }
2047  mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*disVtx.V0vtx) ? mAcc_gfit(*disVtx.V0vtx) : 0;
2048  mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*disVtx.V0vtx) ? mAcc_gmass(*disVtx.V0vtx) : -1;
2049  mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*disVtx.V0vtx) ? mAcc_gmasserr(*disVtx.V0vtx) : -1;
2050  mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*disVtx.V0vtx) ? mAcc_gchisq(*disVtx.V0vtx) : 999999;
2051  mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*disVtx.V0vtx) ? mAcc_gndof(*disVtx.V0vtx) : 0;
2052  mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*disVtx.V0vtx) ? mAcc_gprob(*disVtx.V0vtx) : -1;
2053  trk_px.clear(); trk_py.clear(); trk_pz.clear();
2054  trk_px.push_back( disVtx.p4_V0track1.Px() ); trk_px.push_back( disVtx.p4_V0track2.Px() );
2055  trk_py.push_back( disVtx.p4_V0track1.Py() ); trk_py.push_back( disVtx.p4_V0track2.Py() );
2056  trk_pz.push_back( disVtx.p4_V0track1.Pz() ); trk_pz.push_back( disVtx.p4_V0track2.Pz() );
2057  trk_pxDeco(*cascadeVertices[0]) = trk_px;
2058  trk_pyDeco(*cascadeVertices[0]) = trk_py;
2059  trk_pzDeco(*cascadeVertices[0]) = trk_pz;
2060  trk_px_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Px();
2061  trk_py_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Py();
2062  trk_pz_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Pz();
2063 
2064  result.push_back( fit_result.release() );
2065  }
2066  }
2067  } // m_extraTrk1MassHypo<=0
2068  else { // m_extraTrk1MassHypo>0
2069  std::vector<double> massesJXExtra = massesJX;
2070  massesJXExtra.push_back(m_extraTrk1MassHypo);
2071 
2072  for(const xAOD::TrackParticle* tpExtra : *trackContainer) {
2073  if ( tpExtra->pt()<m_extraTrk1MinPt ) continue;
2074  if ( !m_trkSelector->decision(*tpExtra, nullptr) ) continue;
2075  // Check identical tracks in input
2076  if(std::find(tracksJX.cbegin(),tracksJX.cend(),tpExtra) != tracksJX.cend()) continue;
2077  if(std::find(tracksV0.cbegin(),tracksV0.cend(),tpExtra) != tracksV0.cend()) continue;
2078  if(tpExtra == disVtx.track) continue;
2079 
2080  TLorentzVector tmp;
2081  tmp.SetPtEtaPhiM(tpExtra->pt(),tpExtra->eta(),tpExtra->phi(),m_extraTrk1MassHypo);
2082  if ((p4_moth+tmp).M() < m_MassLower || (p4_moth+tmp).M() > m_MassUpper) continue;
2083 
2084  std::vector<const xAOD::TrackParticle*> tracksJXExtra = tracksJX;
2085  tracksJXExtra.push_back(tpExtra);
2086 
2087  // Apply the user's settings to the fitter
2088  std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState();
2089  // Robustness: http://cdsweb.cern.ch/record/685551
2090  int robustness = 0;
2091  m_iVertexFitter->setRobustness(robustness, *state);
2092  // Build up the topology
2093  // Vertex list
2094  std::vector<Trk::VertexID> vrtList;
2095  std::vector<Trk::VertexID> vrtList2;
2096  // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
2097  // V0 vertex
2098  Trk::VertexID vID1;
2099  if (m_constrV0) {
2100  vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state,m_massV0);
2101  } else {
2102  vID1 = m_iVertexFitter->startVertex(tracksV0,massesV0,*state);
2103  }
2104  vrtList.push_back(vID1);
2105  // Displaced vertex
2106  Trk::VertexID vID2;
2107  if (m_constrDisV) {
2108  vID2 = m_iVertexFitter->nextVertex(tracks3,massesDis3,vrtList,*state,m_massDisV);
2109  } else {
2110  vID2 = m_iVertexFitter->nextVertex(tracks3,massesDis3,vrtList,*state);
2111  }
2112  vrtList2.push_back(vID2);
2113  Trk::VertexID vID3;
2114  if(m_JXSubVtx) {
2115  // JXExtra vertex
2116  vID3 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,*state);
2117  vrtList2.push_back(vID3);
2118  // Mother vertex includes two subvertices (DisV and JX) and extra track
2119  std::vector<const xAOD::TrackParticle*> tp;
2120  std::vector<double> tp_masses;
2121  if(m_constrMainV) {
2122  m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2,*state,m_massMainV);
2123  } else {
2124  m_iVertexFitter->nextVertex(tp,tp_masses,vrtList2,*state);
2125  }
2126  }
2127  else { // m_JXSubVtx=false
2128  // Mother vertex includes just one subvertex (DisV) and JX tracks + extra track
2129  if(m_constrMainV) {
2130  vID3 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList2,*state,m_massMainV);
2131  } else {
2132  vID3 = m_iVertexFitter->nextVertex(tracksJXExtra,massesJXExtra,vrtList2,*state);
2133  }
2134  }
2135  if (m_constrJX && m_jxDaug_num>2) {
2136  std::vector<Trk::VertexID> cnstV;
2137  if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
2138  ATH_MSG_WARNING("addMassConstraint for JX failed");
2139  }
2140  }
2141  if (m_constrJpsi) {
2142  std::vector<Trk::VertexID> cnstV;
2143  if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
2144  ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
2145  }
2146  }
2147  if (m_constrX && m_jxDaug_num==4 && m_massX>0) {
2148  std::vector<Trk::VertexID> cnstV;
2149  if ( !m_iVertexFitter->addMassConstraint(vID3,tracksX,cnstV,*state,m_massX).isSuccess() ) {
2150  ATH_MSG_WARNING("addMassConstraint for X failed");
2151  }
2152  }
2153  // Do the work
2154  std::unique_ptr<Trk::VxCascadeInfo> fit_result = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state) );
2155 
2156  if (fit_result) {
2157  for(auto& v : fit_result->vertices()) {
2158  if(v->nTrackParticles()==0) {
2159  std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
2160  v->setTrackParticleLinks(nullLinkVector);
2161  }
2162  }
2163  // reset links to original tracks
2164  BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
2165 
2166  // necessary to prevent memory leak
2167  fit_result->setSVOwnership(true);
2168 
2169  // Chi2/DOF cut
2170  double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
2171  bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
2172 
2173  const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
2174  const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
2175  size_t iMoth = cascadeVertices.size()-1;
2176  double lxy_SV1_sub = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
2177  double lxy_SV1 = m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[iMoth]);
2178 
2179  if(chi2CutPassed && lxy_SV1>m_lxyDisV_cut && lxy_SV1_sub>m_lxyV0_cut) {
2180  chi2_V0_decor(*cascadeVertices[0]) = disVtx.V0vtx->chiSquared();
2181  ndof_V0_decor(*cascadeVertices[0]) = disVtx.V0vtx->numberDoF();
2182  if(disVtx.V0type==LAMBDA) {
2183  type_V0_decor(*cascadeVertices[0]) = "Lambda";
2184  }
2185  else if(disVtx.V0type==LAMBDABAR) {
2186  type_V0_decor(*cascadeVertices[0]) = "Lambdabar";
2187  }
2188  else if(disVtx.V0type==KS) {
2189  type_V0_decor(*cascadeVertices[0]) = "Ks";
2190  }
2191  mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*disVtx.V0vtx) ? mAcc_gfit(*disVtx.V0vtx) : 0;
2192  mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*disVtx.V0vtx) ? mAcc_gmass(*disVtx.V0vtx) : -1;
2193  mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*disVtx.V0vtx) ? mAcc_gmasserr(*disVtx.V0vtx) : -1;
2194  mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*disVtx.V0vtx) ? mAcc_gchisq(*disVtx.V0vtx) : 999999;
2195  mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*disVtx.V0vtx) ? mAcc_gndof(*disVtx.V0vtx) : 0;
2196  mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*disVtx.V0vtx) ? mAcc_gprob(*disVtx.V0vtx) : -1;
2197  trk_px.clear(); trk_py.clear(); trk_pz.clear();
2198  trk_px.push_back( disVtx.p4_V0track1.Px() ); trk_px.push_back( disVtx.p4_V0track2.Px() );
2199  trk_py.push_back( disVtx.p4_V0track1.Py() ); trk_py.push_back( disVtx.p4_V0track2.Py() );
2200  trk_pz.push_back( disVtx.p4_V0track1.Pz() ); trk_pz.push_back( disVtx.p4_V0track2.Pz() );
2201  trk_pxDeco(*cascadeVertices[0]) = trk_px;
2202  trk_pyDeco(*cascadeVertices[0]) = trk_py;
2203  trk_pzDeco(*cascadeVertices[0]) = trk_pz;
2204  trk_px_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Px();
2205  trk_py_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Py();
2206  trk_pz_deco(*cascadeVertices[1]) = disVtx.p4_disVtrack.Pz();
2207 
2208  result.push_back( fit_result.release() );
2209  }
2210  }
2211  } // loop over trackContainer
2212  } // m_extraTrk1MassHypo>0
2213 
2214  return result;
2215  }
2216 
2217  void JpsiXPlusDisplaced::fitV0Container(xAOD::VertexContainer* V0ContainerNew, const std::vector<const xAOD::TrackParticle*>& selectedTracks, const std::vector<const xAOD::TrackParticleContainer*>& trackCols) const {
2218  const EventContext& ctx = Gaudi::Hive::currentContext();
2219 
2220  SG::AuxElement::Decorator<std::string> mDec_type("Type_V0Vtx");
2221  SG::AuxElement::Decorator<int> mDec_gfit("gamma_fit");
2222  SG::AuxElement::Decorator<float> mDec_gmass("gamma_mass");
2223  SG::AuxElement::Decorator<float> mDec_gmasserr("gamma_massError");
2224  SG::AuxElement::Decorator<float> mDec_gchisq("gamma_chisq");
2225  SG::AuxElement::Decorator<int> mDec_gndof("gamma_ndof");
2226  SG::AuxElement::Decorator<float> mDec_gprob("gamma_probability");
2227 
2228  std::vector<const xAOD::TrackParticle*> posTracks;
2229  std::vector<const xAOD::TrackParticle*> negTracks;
2230  for(const xAOD::TrackParticle* TP : selectedTracks) {
2231  if(TP->charge()>0) posTracks.push_back(TP);
2232  else negTracks.push_back(TP);
2233  }
2234 
2235  for(const xAOD::TrackParticle* TP1 : posTracks) {
2236  const Trk::Perigee& aPerigee1 = TP1->perigeeParameters();
2237  for(const xAOD::TrackParticle* TP2 : negTracks) {
2238  const Trk::Perigee& aPerigee2 = TP2->perigeeParameters();
2239  int sflag(0), errorcode(0);
2240  Amg::Vector3D startingPoint = m_vertexEstimator->getCirclesIntersectionPoint(&aPerigee1,&aPerigee2,sflag,errorcode);
2241  if (errorcode != 0) {startingPoint(0) = 0.0; startingPoint(1) = 0.0; startingPoint(2) = 0.0;}
2242 
2243  if (errorcode == 0 || errorcode == 5 || errorcode == 6 || errorcode == 8) {
2244  Trk::PerigeeSurface perigeeSurface(startingPoint);
2245  const Trk::TrackParameters* extrapolatedPerigee1 = m_extrapolator->extrapolate(ctx,TP1->perigeeParameters(), perigeeSurface).release();
2246  const Trk::TrackParameters* extrapolatedPerigee2 = m_extrapolator->extrapolate(ctx,TP2->perigeeParameters(), perigeeSurface).release();
2247  std::vector<std::unique_ptr<const Trk::TrackParameters> > cleanup;
2248  if(!extrapolatedPerigee1) extrapolatedPerigee1 = &TP1->perigeeParameters();
2249  else cleanup.push_back(std::unique_ptr<const Trk::TrackParameters>(extrapolatedPerigee1));
2250  if(!extrapolatedPerigee2) extrapolatedPerigee2 = &TP2->perigeeParameters();
2251  else cleanup.push_back(std::unique_ptr<const Trk::TrackParameters>(extrapolatedPerigee2));
2252  if(extrapolatedPerigee1 && extrapolatedPerigee2) {
2253  bool pass = false;
2254  TLorentzVector v1; TLorentzVector v2;
2255  if(!pass) {
2256  v1.SetXYZM(extrapolatedPerigee1->momentum().x(),extrapolatedPerigee1->momentum().y(),extrapolatedPerigee1->momentum().z(),m_mass_proton);
2257  v2.SetXYZM(extrapolatedPerigee2->momentum().x(),extrapolatedPerigee2->momentum().y(),extrapolatedPerigee2->momentum().z(),m_mass_pion);
2258  if((v1+v2).M()>900.0 && (v1+v2).M()<1350.0) pass = true;
2259  }
2260  if(!pass) {
2261  v1.SetXYZM(extrapolatedPerigee1->momentum().x(),extrapolatedPerigee1->momentum().y(),extrapolatedPerigee1->momentum().z(),m_mass_pion);
2262  v2.SetXYZM(extrapolatedPerigee2->momentum().x(),extrapolatedPerigee2->momentum().y(),extrapolatedPerigee2->momentum().z(),m_mass_proton);
2263  if((v1+v2).M()>900.0 && (v1+v2).M()<1350.0) pass = true;
2264  }
2265  if(!pass) {
2266  v1.SetXYZM(extrapolatedPerigee1->momentum().x(),extrapolatedPerigee1->momentum().y(),extrapolatedPerigee1->momentum().z(),m_mass_pion);
2267  v2.SetXYZM(extrapolatedPerigee2->momentum().x(),extrapolatedPerigee2->momentum().y(),extrapolatedPerigee2->momentum().z(),m_mass_pion);
2268  if((v1+v2).M()>300.0 && (v1+v2).M()<700.0) pass = true;
2269  }
2270  if(pass) {
2271  std::vector<const xAOD::TrackParticle*> tracksV0;
2272  tracksV0.push_back(TP1); tracksV0.push_back(TP2);
2273  std::unique_ptr<xAOD::Vertex> V0vtx = std::unique_ptr<xAOD::Vertex>( m_iV0Fitter->fit(tracksV0, startingPoint) );
2274  if(V0vtx && V0vtx->chiSquared()>=0) {
2275  double chi2DOF = V0vtx->chiSquared()/V0vtx->numberDoF();
2276  if(chi2DOF>m_chi2cut_V0) continue;
2277 
2278  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);
2279  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);
2280  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);
2281  if(massSig_V0_Lambda1<=massSig_V0_Lambda2 && massSig_V0_Lambda1<=massSig_V0_Ks) {
2282  mDec_type(*V0vtx.get()) = "Lambda";
2283  }
2284  else if(massSig_V0_Lambda2<=massSig_V0_Lambda1 && massSig_V0_Lambda2<=massSig_V0_Ks) {
2285  mDec_type(*V0vtx.get()) = "Lambdabar";
2286  }
2287  else if(massSig_V0_Ks<=massSig_V0_Lambda1 && massSig_V0_Ks<=massSig_V0_Lambda2) {
2288  mDec_type(*V0vtx.get()) = "Ks";
2289  }
2290 
2291  int gamma_fit = 0; int gamma_ndof = 0; double gamma_chisq = 999999.;
2292  double gamma_prob = -1., gamma_mass = -1., gamma_massErr = -1.;
2293  std::unique_ptr<xAOD::Vertex> gammaVtx = std::unique_ptr<xAOD::Vertex>( m_iGammaFitter->fit(tracksV0, m_V0Tools->vtx(V0vtx.get())) );
2294  if (gammaVtx) {
2295  gamma_fit = 1;
2296  gamma_mass = m_V0Tools->invariantMass(gammaVtx.get(),m_mass_e,m_mass_e);
2297  gamma_massErr = m_V0Tools->invariantMassError(gammaVtx.get(),m_mass_e,m_mass_e);
2298  gamma_chisq = m_V0Tools->chisq(gammaVtx.get());
2299  gamma_ndof = m_V0Tools->ndof(gammaVtx.get());
2300  gamma_prob = m_V0Tools->vertexProbability(gammaVtx.get());
2301  }
2302  mDec_gfit(*V0vtx.get()) = gamma_fit;
2303  mDec_gmass(*V0vtx.get()) = gamma_mass;
2304  mDec_gmasserr(*V0vtx.get()) = gamma_massErr;
2305  mDec_gchisq(*V0vtx.get()) = gamma_chisq;
2306  mDec_gndof(*V0vtx.get()) = gamma_ndof;
2307  mDec_gprob(*V0vtx.get()) = gamma_prob;
2308 
2309  xAOD::BPhysHelper V0_helper(V0vtx.get());
2310  V0_helper.setRefTrks(); // AOD only method
2311 
2312  if(not trackCols.empty()){
2313  try {
2314  JpsiUpsilonCommon::RelinkVertexTracks(trackCols, V0vtx.get());
2315  } catch (std::runtime_error const& e) {
2316  ATH_MSG_ERROR(e.what());
2317  return;
2318  }
2319  }
2320 
2321  V0ContainerNew->push_back(std::move(V0vtx));
2322  }
2323  }
2324  }
2325  }
2326  }
2327  }
2328  }
2329 
2330  template<size_t NTracks>
2331  const xAOD::Vertex* JpsiXPlusDisplaced::FindVertex(const xAOD::VertexContainer* cont, const xAOD::Vertex* v) const {
2332  for (const xAOD::Vertex* v1 : *cont) {
2333  assert(v1->nTrackParticles() == NTracks);
2334  std::array<const xAOD::TrackParticle*, NTracks> a1;
2335  std::array<const xAOD::TrackParticle*, NTracks> a2;
2336  for(size_t i=0; i<NTracks; i++){
2337  a1[i] = v1->trackParticle(i);
2338  a2[i] = v->trackParticle(i);
2339  }
2340  std::sort(a1.begin(), a1.end());
2341  std::sort(a2.begin(), a2.end());
2342  if(a1 == a2) return v1;
2343  }
2344  return nullptr;
2345  }
2346 }
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:848
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:794
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:228
Trk::d0
@ d0
Definition: ParamDefs.h:63
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
DerivationFramework::BPhysPVCascadeTools
Definition: BPhysPVCascadeTools.h:34
DerivationFramework::JpsiXPlusDisplaced::MesonCandidate
Definition: JpsiXPlusDisplaced.h: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 Sun Dec 22 2024 21:12:48 for ATLAS Offline Software by   doxygen 1.8.18