ATLAS Offline Software
Loading...
Searching...
No Matches
JpsiXPlus2V0.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3 Contact: Xin Chen <xin.chen@cern.ch>
4*/
5#include "JpsiXPlus2V0.h"
12#include "BPhysPVCascadeTools.h"
17#include "VxVertex/RecVertex.h"
20#include <algorithm>
21#include <functional>
22
23namespace DerivationFramework {
24 typedef ElementLink<xAOD::VertexContainer> VertexLink;
25 typedef std::vector<VertexLink> VertexLinkVector;
26
27 using Analysis::JpsiUpsilonCommon;
28
29 JpsiXPlus2V0::JpsiXPlus2V0(const std::string& type, const std::string& name, const IInterface* parent) : base_class(type,name,parent),
30 m_vertexJXContainerKey("InputJXVertices"),
32 m_cascadeOutputKeys({"JpsiXPlus2V0_SubVtx1", "JpsiXPlus2V0_SubVtx2", "JpsiXPlus2V0_SubVtx3", "JpsiXPlus2V0_MainVtx"}),
33 m_v0VtxOutputKey(""),
34 m_TrkParticleCollection("InDetTrackParticles"),
35 m_VxPrimaryCandidateName("PrimaryVertices"),
36 m_pvContainerName("PrimaryVertices"),
37 m_refPVContainerName("RefittedPrimaryVertices"),
38 m_eventInfo_key("EventInfo"),
39 m_RelinkContainers({"InDetTrackParticles","InDetLargeD0TrackParticles"}),
40 m_useImprovedMass(false),
41 m_jxMassLower(0.0),
42 m_jxMassUpper(30000.0),
43 m_jpsiMassLower(0.0),
44 m_jpsiMassUpper(20000.0),
45 m_diTrackMassLower(-1.0),
46 m_diTrackMassUpper(-1.0),
47 m_V01Hypothesis("Ks"),
48 m_V02Hypothesis("Lambda"),
49 m_LambdaMassLower(0.0),
50 m_LambdaMassUpper(10000.0),
51 m_KsMassLower(0.0),
52 m_KsMassUpper(10000.0),
53 m_lxyV01_cut(-999.0),
54 m_lxyV02_cut(-999.0),
55 m_minMass_gamma(-1.0),
56 m_chi2cut_gamma(-1.0),
57 m_JXV02MassLower(0.0),
58 m_JXV02MassUpper(30000.0),
59 m_MassLower(0.0),
60 m_MassUpper(31000.0),
61 m_jxDaug_num(4),
62 m_jxDaug1MassHypo(-1),
63 m_jxDaug2MassHypo(-1),
64 m_jxDaug3MassHypo(-1),
65 m_jxDaug4MassHypo(-1),
66 m_massJX(-1),
67 m_massJpsi(-1),
68 m_massX(-1),
69 m_massLd(-1),
70 m_massKs(-1),
71 m_massJXV02(-1),
72 m_massMainV(-1),
73 m_constrJX(false),
74 m_constrJpsi(false),
75 m_constrX(false),
76 m_constrV01(false),
77 m_constrV02(false),
78 m_constrJXV02(false),
79 m_constrMainV(false),
80 m_cascadeFitWithPV(0),
81 m_firstDecayAtPV(false),
82 m_JXSubVtx(false),
83 m_JXV02SubVtx(false),
84 m_chi2cut_JX(-1.0),
85 m_chi2cut_V0(-1.0),
86 m_chi2cut(-1.0),
87 m_useTRT(false),
88 m_ptTRT(450),
89 m_d0_cut(2),
90 m_maxJXCandidates(0),
91 m_maxV0Candidates(0),
92 m_maxMainVCandidates(0),
93 m_iVertexFitter("Trk::TrkVKalVrtFitter"),
94 m_iV0Fitter("Trk::V0VertexFitter"),
95 m_iGammaFitter("Trk::TrkVKalVrtFitter"),
96 m_pvRefitter("Analysis::PrimaryVertexRefitter"),
97 m_V0Tools("Trk::V0Tools"),
98 m_trackToVertexTool("Reco::TrackToVertex"),
99 m_v0TrkSelector("InDet::TrackSelectorTool"),
100 m_CascadeTools("DerivationFramework::CascadeTools"),
101 m_vertexEstimator("InDet::VertexPointEstimator"),
102 m_extrapolator("Trk::Extrapolator/AtlasExtrapolator")
103 {
104 declareProperty("JXVertices", m_vertexJXContainerKey);
105 declareProperty("V0Vertices", m_vertexV0ContainerKey);
106 declareProperty("JXVtxHypoNames", m_vertexJXHypoNames);
107 declareProperty("CascadeVertexCollections", m_cascadeOutputKeys); // size is 3 or 4 only
108 declareProperty("OutoutV0VtxCollection", m_v0VtxOutputKey);
109 declareProperty("TrackParticleCollection", m_TrkParticleCollection);
110 declareProperty("VxPrimaryCandidateName", m_VxPrimaryCandidateName);
111 declareProperty("PVContainerName", m_pvContainerName);
112 declareProperty("RefPVContainerName", m_refPVContainerName);
113 declareProperty("EventInfoKey", m_eventInfo_key);
114 declareProperty("RelinkTracks", m_RelinkContainers);
115 declareProperty("UseImprovedMass", m_useImprovedMass);
116 declareProperty("JXMassLowerCut", m_jxMassLower); // only effective when m_jxDaug_num>2
117 declareProperty("JXMassUpperCut", m_jxMassUpper); // only effective when m_jxDaug_num>2
118 declareProperty("JpsiMassLowerCut", m_jpsiMassLower);
119 declareProperty("JpsiMassUpperCut", m_jpsiMassUpper);
120 declareProperty("DiTrackMassLower", m_diTrackMassLower); // only effective when m_jxDaug_num=4
121 declareProperty("DiTrackMassUpper", m_diTrackMassUpper); // only effective when m_jxDaug_num=4
122 declareProperty("V01Hypothesis", m_V01Hypothesis); // "Ks", "Lambda" or "Lambda/Ks"
123 declareProperty("V02Hypothesis", m_V02Hypothesis); // "Ks", "Lambda" or "Lambda/Ks"
124 declareProperty("LxyV01Cut", m_lxyV01_cut);
125 declareProperty("LxyV02Cut", m_lxyV02_cut);
126 declareProperty("LambdaMassLowerCut", m_LambdaMassLower);
127 declareProperty("LambdaMassUpperCut", m_LambdaMassUpper);
128 declareProperty("KsMassLowerCut", m_KsMassLower);
129 declareProperty("KsMassUpperCut", m_KsMassUpper);
130 declareProperty("MassCutGamma", m_minMass_gamma);
131 declareProperty("Chi2CutGamma", m_chi2cut_gamma);
132 declareProperty("JXV02MassLowerCut", m_JXV02MassLower); // only effective when m_JXSubVtx=true & m_JXV02SubVtx=true
133 declareProperty("JXV02MassUpperCut", m_JXV02MassUpper); // only effective when m_JXSubVtx=true & m_JXV02SubVtx=true
134 declareProperty("MassLowerCut", m_MassLower);
135 declareProperty("MassUpperCut", m_MassUpper);
136 declareProperty("HypothesisName", m_hypoName = "TQ");
137 declareProperty("NumberOfJXDaughters", m_jxDaug_num); // 2, or 3, or 4 only
138 declareProperty("JXDaug1MassHypo", m_jxDaug1MassHypo);
139 declareProperty("JXDaug2MassHypo", m_jxDaug2MassHypo);
140 declareProperty("JXDaug3MassHypo", m_jxDaug3MassHypo);
141 declareProperty("JXDaug4MassHypo", m_jxDaug4MassHypo);
142 declareProperty("JXMass", m_massJX); // only effective when m_jxDaug_num>2
143 declareProperty("JpsiMass", m_massJpsi);
144 declareProperty("XMass", m_massX); // only effective when m_jxDaug_num=4
145 declareProperty("LambdaMass", m_massLd);
146 declareProperty("KsMass", m_massKs);
147 declareProperty("JXV02VtxMass", m_massJXV02); // mass of JX + 2nd V0
148 declareProperty("MainVtxMass", m_massMainV);
149 declareProperty("ApplyJXMassConstraint", m_constrJX); // only effective when m_jxDaug_num>2
150 declareProperty("ApplyJpsiMassConstraint", m_constrJpsi);
151 declareProperty("ApplyXMassConstraint", m_constrX); // only effective when m_jxDaug_num=4
152 declareProperty("ApplyV01MassConstraint", m_constrV01); // first V0
153 declareProperty("ApplyV02MassConstraint", m_constrV02); // second V0
154 declareProperty("ApplyJXV02MassConstraint", m_constrJXV02); // constrain JX + 2nd V0
155 declareProperty("ApplyMainVMassConstraint", m_constrMainV);
156 declareProperty("DoCascadeFitWithPV", m_cascadeFitWithPV);
157 declareProperty("FirstDecayAtPV", m_firstDecayAtPV);
158 declareProperty("HasJXSubVertex", m_JXSubVtx);
159 declareProperty("HasJXV02SubVertex", m_JXV02SubVtx); // only effective when m_JXSubVtx=true
160 declareProperty("Chi2CutJX", m_chi2cut_JX);
161 declareProperty("Chi2CutV0", m_chi2cut_V0);
162 declareProperty("Chi2Cut", m_chi2cut);
163 declareProperty("UseTRT", m_useTRT);
164 declareProperty("PtTRT", m_ptTRT);
165 declareProperty("Trackd0Cut", m_d0_cut);
166 declareProperty("MaxJXCandidates", m_maxJXCandidates);
167 declareProperty("MaxV0Candidates", m_maxV0Candidates);
168 declareProperty("MaxMainVCandidates", m_maxMainVCandidates);
169 declareProperty("RefitPV", m_refitPV = true);
170 declareProperty("MaxnPV", m_PV_max = 1000);
171 declareProperty("MinNTracksInPV", m_PV_minNTracks = 0);
172 declareProperty("DoVertexType", m_DoVertexType = 7);
173 declareProperty("TrkVertexFitterTool", m_iVertexFitter);
174 declareProperty("V0VertexFitterTool", m_iV0Fitter);
175 declareProperty("GammaFitterTool", m_iGammaFitter);
176 declareProperty("PVRefitter", m_pvRefitter);
177 declareProperty("V0Tools", m_V0Tools);
178 declareProperty("TrackToVertexTool", m_trackToVertexTool);
179 declareProperty("V0TrackSelectorTool", m_v0TrkSelector);
180 declareProperty("CascadeTools", m_CascadeTools);
181 declareProperty("VertexPointEstimator", m_vertexEstimator);
182 declareProperty("Extrapolator", m_extrapolator);
183 }
184
186 if((m_V01Hypothesis != "Ks" && m_V01Hypothesis != "Lambda" && m_V01Hypothesis != "Lambda/Ks" && m_V01Hypothesis != "Ks/Lambda") ||
187 (m_V02Hypothesis != "Ks" && m_V02Hypothesis != "Lambda" && m_V02Hypothesis != "Lambda/Ks" && m_V02Hypothesis != "Ks/Lambda")) {
188 ATH_MSG_FATAL("Incorrect V0 container hypothesis - not recognized");
189 return StatusCode::FAILURE;
190 }
191
193 ATH_MSG_FATAL("Incorrect number of JX daughters");
194 return StatusCode::FAILURE;
195 }
196
197 if(m_vertexV0ContainerKey.key()=="" && m_v0VtxOutputKey.key()=="") {
198 ATH_MSG_FATAL("Input and output V0 container names can not be both empty");
199 return StatusCode::FAILURE;
200 }
201
202 // retrieving vertex Fitter
203 ATH_CHECK( m_iVertexFitter.retrieve() );
204
205 // retrieving V0 vertex Fitter
206 ATH_CHECK( m_iV0Fitter.retrieve() );
207
208 // retrieving photon conversion vertex Fitter
209 ATH_CHECK( m_iGammaFitter.retrieve() );
210
211 // retrieving primary vertex refitter
212 ATH_CHECK( m_pvRefitter.retrieve() );
213
214 // retrieving the V0 tool
215 ATH_CHECK( m_V0Tools.retrieve() );
216
217 // retrieving the TrackToVertex extrapolator tool
218 ATH_CHECK( m_trackToVertexTool.retrieve() );
219
220 // retrieving the V0 track selector tool
221 ATH_CHECK( m_v0TrkSelector.retrieve() );
222
223 // retrieving the Cascade tools
224 ATH_CHECK( m_CascadeTools.retrieve() );
225
226 // retrieving the vertex point estimator
227 ATH_CHECK( m_vertexEstimator.retrieve() );
228
229 // retrieving the extrapolator
230 ATH_CHECK( m_extrapolator.retrieve() );
231
232 ATH_CHECK( m_vertexJXContainerKey.initialize() );
234 ATH_CHECK( m_VxPrimaryCandidateName.initialize() );
235 ATH_CHECK( m_TrkParticleCollection.initialize() );
236 ATH_CHECK( m_pvContainerName.initialize() );
237 ATH_CHECK( m_refPVContainerName.initialize() );
238 ATH_CHECK( m_cascadeOutputKeys.initialize() );
239 ATH_CHECK( m_eventInfo_key.initialize() );
240 ATH_CHECK( m_RelinkContainers.initialize() );
242
243 auto gendata = std::make_shared<GenData>();
244
245 m_mass_e = gendata->particleMass(MC::ELECTRON).value();
246 m_mass_mu = gendata->particleMass(MC::MUON).value();
247 m_mass_pion = gendata->particleMass(MC::PIPLUS).value();
248 m_mass_proton = gendata->particleMass(MC::PROTON).value();
249 m_mass_Lambda = gendata->particleMass(MC::LAMBDA0).value();
250 m_mass_Lambdab = gendata->particleMass(MC::LAMBDAB0).value();
251 m_mass_BCPLUS = gendata->particleMass(MC::BCPLUS).value();
252 m_mass_Ks = gendata->particleMass(MC::K0S).value();
253 m_mass_Bpm = gendata->particleMass(521).value();
254 m_mass_phi = gendata->particleMass(333).value();
255
256 m_massesV0_ppi.push_back(m_mass_proton);
257 m_massesV0_ppi.push_back(m_mass_pion);
258 m_massesV0_pip.push_back(m_mass_pion);
259 m_massesV0_pip.push_back(m_mass_proton);
260 m_massesV0_pipi.push_back(m_mass_pion);
261 m_massesV0_pipi.push_back(m_mass_pion);
262
263 // retrieve particle masses
264 if(m_constrJpsi && m_massJpsi<0) m_massJpsi = gendata->particleMass(MC::JPSI).value();
265 if(m_jxDaug_num>=3 && m_constrJX && m_massJX<0) m_massJX = gendata->particleMass(MC::PSI2S).value();
267 if(m_constrV01 || m_constrV02) {
270 }
273
278
279 return StatusCode::SUCCESS;
280 }
281
282 StatusCode JpsiXPlus2V0::performSearch(std::vector<Trk::VxCascadeInfo*>& cascadeinfoContainer, const std::vector<std::pair<const xAOD::Vertex*,V0Enum> >& selectedV0Candidates, const EventContext& ctx) const {
283 ATH_MSG_DEBUG( "JpsiXPlus2V0::performSearch" );
284 if(selectedV0Candidates.size()==0) return StatusCode::SUCCESS;
285
286 // Get all track containers when m_RelinkContainers is not empty
287 std::vector<const xAOD::TrackParticleContainer*> trackCols;
290 ATH_CHECK( handle.isValid() );
291 trackCols.push_back(handle.cptr());
292 }
293
294 // Get default PV container
296 ATH_CHECK( defaultPVContainer.isValid() );
297
298 // Get PV container
300 ATH_CHECK( pvContainer.isValid() );
301
302 // Get Jpsi+X container
304 ATH_CHECK( jxContainer.isValid() );
305
306 std::vector<double> massesJX{m_jxDaug1MassHypo, m_jxDaug2MassHypo};
307 if(m_jxDaug_num>=3) massesJX.push_back(m_jxDaug3MassHypo);
308 if(m_jxDaug_num==4) massesJX.push_back(m_jxDaug4MassHypo);
309
310 // Select the JX candidates before calling cascade fit
311 std::vector<const xAOD::Vertex*> selectedJXCandidates;
312 for(auto vxcItr=jxContainer.ptr()->begin(); vxcItr!=jxContainer.ptr()->end(); ++vxcItr) {
313 // Check the passed flag first
314 const xAOD::Vertex* vtx = *vxcItr;
315 bool passed = false;
316 for(const std::string& name : m_vertexJXHypoNames) {
317 SG::Accessor<Char_t> flagAcc("passed_"+name);
318 if(flagAcc.isAvailable(*vtx) && flagAcc(*vtx)) {
319 passed = true;
320 }
321 }
322 if(m_vertexJXHypoNames.size() && !passed) continue;
323
324 // Add loose cut on Jpsi mass from e.g. JX -> Jpsi pi+ pi-
325 TLorentzVector p4_mu1, p4_mu2;
326 p4_mu1.SetPtEtaPhiM(vtx->trackParticle(0)->pt(),vtx->trackParticle(0)->eta(),vtx->trackParticle(0)->phi(), m_jxDaug1MassHypo);
327 p4_mu2.SetPtEtaPhiM(vtx->trackParticle(1)->pt(),vtx->trackParticle(1)->eta(),vtx->trackParticle(1)->phi(), m_jxDaug2MassHypo);
328 double mass_jpsi = (p4_mu1 + p4_mu2).M();
329 if (mass_jpsi < m_jpsiMassLower || mass_jpsi > m_jpsiMassUpper) continue;
330
331 TLorentzVector p4_trk1, p4_trk2;
332 if(m_jxDaug_num>=3) p4_trk1.SetPtEtaPhiM(vtx->trackParticle(2)->pt(),vtx->trackParticle(2)->eta(),vtx->trackParticle(2)->phi(), m_jxDaug3MassHypo);
333 if(m_jxDaug_num==4) p4_trk2.SetPtEtaPhiM(vtx->trackParticle(3)->pt(),vtx->trackParticle(3)->eta(),vtx->trackParticle(3)->phi(), m_jxDaug4MassHypo);
334
335 if(m_jxDaug_num==3) {
336 double mass_jx = (p4_mu1 + p4_mu2 + p4_trk1).M();
337 if(m_useImprovedMass && m_massJpsi>0) mass_jx += - (p4_mu1 + p4_mu2).M() + m_massJpsi;
338 if(mass_jx < m_jxMassLower || mass_jx > m_jxMassUpper) continue;
339 }
340 else if(m_jxDaug_num==4) {
341 double mass_jx = (p4_mu1 + p4_mu2 + p4_trk1 + p4_trk2).M();
342 if(m_useImprovedMass && m_massJpsi>0) mass_jx += - (p4_mu1 + p4_mu2).M() + m_massJpsi;
343 if(mass_jx < m_jxMassLower || mass_jx > m_jxMassUpper) continue;
344
346 double mass_diTrk = (p4_trk1 + p4_trk2).M();
347 if(mass_diTrk < m_diTrackMassLower || mass_diTrk > m_diTrackMassUpper) continue;
348 }
349 }
350
351 double chi2DOF = vtx->chiSquared()/vtx->numberDoF();
352 if(m_chi2cut_JX>0 && chi2DOF>m_chi2cut_JX) continue;
353
354 selectedJXCandidates.push_back(vtx);
355 }
356 if(selectedJXCandidates.size()==0) return StatusCode::SUCCESS;
357
358 std::sort( selectedJXCandidates.begin(), selectedJXCandidates.end(), [](const xAOD::Vertex* a, const xAOD::Vertex* b) { return a->chiSquared()/a->numberDoF() < b->chiSquared()/b->numberDoF(); } );
359 if(m_maxJXCandidates>0 && selectedJXCandidates.size()>m_maxJXCandidates) {
360 selectedJXCandidates.erase(selectedJXCandidates.begin()+m_maxJXCandidates, selectedJXCandidates.end());
361 }
362
363 // Select JX+V0+V0 candidates
364 std::vector<const xAOD::TrackParticle*> tracksJX; tracksJX.reserve(m_jxDaug_num);
365 std::vector<const xAOD::TrackParticle*> tracksV01; tracksV01.reserve(2);
366 for(auto jxItr=selectedJXCandidates.cbegin(); jxItr!=selectedJXCandidates.cend(); ++jxItr) {
367 tracksJX.clear();
368 for(size_t i=0; i<(*jxItr)->nTrackParticles(); i++) tracksJX.push_back((*jxItr)->trackParticle(i));
369 for(auto V0Itr1=selectedV0Candidates.cbegin(); V0Itr1!=selectedV0Candidates.cend(); ++V0Itr1) {
370 if(std::find(tracksJX.cbegin(), tracksJX.cend(), V0Itr1->first->trackParticle(0)) != tracksJX.cend()) continue;
371 if(std::find(tracksJX.cbegin(), tracksJX.cend(), V0Itr1->first->trackParticle(1)) != tracksJX.cend()) continue;
372 tracksV01.clear();
373 for(size_t j=0; j<V0Itr1->first->nTrackParticles(); j++) tracksV01.push_back(V0Itr1->first->trackParticle(j));
374 for(auto V0Itr2=V0Itr1+1; V0Itr2!=selectedV0Candidates.cend(); ++V0Itr2) {
375 if(std::find(tracksJX.cbegin(), tracksJX.cend(), V0Itr2->first->trackParticle(0)) != tracksJX.cend()) continue;
376 if(std::find(tracksJX.cbegin(), tracksJX.cend(), V0Itr2->first->trackParticle(1)) != tracksJX.cend()) continue;
377 if(std::find(tracksV01.cbegin(), tracksV01.cend(), V0Itr2->first->trackParticle(0)) != tracksV01.cend()) continue;
378 if(std::find(tracksV01.cbegin(), tracksV01.cend(), V0Itr2->first->trackParticle(1)) != tracksV01.cend()) continue;
379
380 int numberOfVertices = 0;
381 if((m_V01Hypothesis=="Lambda/Ks" || m_V01Hypothesis=="Ks/Lambda") && (m_V02Hypothesis=="Lambda/Ks" || m_V02Hypothesis=="Ks/Lambda")) {
382 numberOfVertices = m_JXSubVtx && m_JXV02SubVtx ? 2 : 1;
383 }
384 else if((m_V01Hypothesis == "Ks" || m_V01Hypothesis == "Lambda") && (m_V02Hypothesis=="Lambda/Ks" || m_V02Hypothesis=="Ks/Lambda")) {
385 if((m_V01Hypothesis == "Ks" && V0Itr1->second==KS) ||
386 (m_V01Hypothesis == "Lambda" && (V0Itr1->second==LAMBDA || V0Itr1->second==LAMBDABAR))) numberOfVertices = 1;
387 if((m_V01Hypothesis == "Ks" && V0Itr2->second==KS) ||
388 (m_V01Hypothesis == "Lambda" && (V0Itr2->second==LAMBDA || V0Itr2->second==LAMBDABAR))) {
389 if(numberOfVertices == 1) numberOfVertices = m_JXSubVtx && m_JXV02SubVtx ? 2 : 1;
390 else numberOfVertices = -1;
391 }
392 }
393 else if((m_V02Hypothesis == "Ks" || m_V02Hypothesis == "Lambda") && (m_V01Hypothesis=="Lambda/Ks" || m_V01Hypothesis=="Ks/Lambda")) {
394 if((m_V02Hypothesis == "Ks" && V0Itr2->second==KS) ||
395 (m_V02Hypothesis == "Lambda" && (V0Itr2->second==LAMBDA || V0Itr2->second==LAMBDABAR))) numberOfVertices = 1;
396 if((m_V02Hypothesis == "Ks" && V0Itr1->second==KS) ||
397 (m_V02Hypothesis == "Lambda" && (V0Itr1->second==LAMBDA || V0Itr1->second==LAMBDABAR))) {
398 if(numberOfVertices == 1) numberOfVertices = m_JXSubVtx && m_JXV02SubVtx ? 2 : 1;
399 else numberOfVertices = -1;
400 }
401 }
402 else if(m_V01Hypothesis == "Ks" && m_V02Hypothesis == "Lambda") {
403 if(V0Itr1->second==KS && (V0Itr2->second==LAMBDA || V0Itr2->second==LAMBDABAR)) numberOfVertices = 1;
404 else if(V0Itr2->second==KS && (V0Itr1->second==LAMBDA || V0Itr1->second==LAMBDABAR)) numberOfVertices = -1;
405 }
406 else if(m_V01Hypothesis == "Lambda" && m_V02Hypothesis == "Ks") {
407 if((V0Itr1->second==LAMBDA || V0Itr1->second==LAMBDABAR) && V0Itr2->second==KS) numberOfVertices = 1;
408 else if((V0Itr2->second==LAMBDA || V0Itr2->second==LAMBDABAR) && V0Itr1->second==KS) numberOfVertices = -1;
409 }
410 else if(m_V01Hypothesis == "Ks" && m_V02Hypothesis == "Ks") {
411 if(V0Itr1->second==KS && V0Itr2->second==KS) numberOfVertices = m_JXSubVtx && m_JXV02SubVtx ? 2 : 1;
412 }
413 else if(m_V01Hypothesis == "Lambda" && m_V02Hypothesis == "Lambda") {
414 if((V0Itr1->second==LAMBDA || V0Itr1->second==LAMBDABAR) && (V0Itr2->second==LAMBDA || V0Itr2->second==LAMBDABAR)) numberOfVertices = m_JXSubVtx && m_JXV02SubVtx ? 2 : 1;
415 }
416
417 if(numberOfVertices==2) {
418 Trk::VxCascadeInfo* result1 = fitMainVtx(ctx, *jxItr, massesJX, V0Itr1->first, V0Itr1->second, V0Itr2->first, V0Itr2->second, trackCols, defaultPVContainer.cptr(), pvContainer.cptr());
419 if(result1) cascadeinfoContainer.push_back(result1);
420 Trk::VxCascadeInfo* result2 = fitMainVtx(ctx, *jxItr, massesJX, V0Itr2->first, V0Itr2->second, V0Itr1->first, V0Itr1->second, trackCols, defaultPVContainer.cptr(), pvContainer.cptr());
421 if(result2) cascadeinfoContainer.push_back(result2);
422 }
423 else if(numberOfVertices==1) {
424 Trk::VxCascadeInfo* result = fitMainVtx(ctx, *jxItr, massesJX, V0Itr1->first, V0Itr1->second, V0Itr2->first, V0Itr2->second, trackCols, defaultPVContainer.cptr(), pvContainer.cptr());
425 if(result) cascadeinfoContainer.push_back(result);
426 }
427 else if(numberOfVertices==-1) {
428 Trk::VxCascadeInfo* result = fitMainVtx(ctx, *jxItr, massesJX, V0Itr2->first, V0Itr2->second, V0Itr1->first, V0Itr1->second, trackCols, defaultPVContainer.cptr(), pvContainer.cptr());
429 if(result) cascadeinfoContainer.push_back(result);
430 }
431 }
432 }
433 }
434
435 return StatusCode::SUCCESS;
436 }
437
438 StatusCode JpsiXPlus2V0::addBranches(const EventContext& ctx) const {
439 size_t topoN = 4;
440 if(!m_JXSubVtx) topoN--;
441
442 if(m_cascadeOutputKeys.size() != topoN) {
443 ATH_MSG_FATAL("Incorrect number of output cascade vertices");
444 return StatusCode::FAILURE;
445 }
446
447 std::array<SG::WriteHandle<xAOD::VertexContainer>, 4> VtxWriteHandles; int ikey(0);
449 VtxWriteHandles[ikey] = SG::WriteHandle<xAOD::VertexContainer>(key, ctx);
450 ATH_CHECK( VtxWriteHandles[ikey].record(std::make_unique<xAOD::VertexContainer>(), std::make_unique<xAOD::VertexAuxContainer>()) );
451 ikey++;
452 }
453
454 //----------------------------------------------------
455 // retrieve primary vertices
456 //----------------------------------------------------
457 const xAOD::Vertex* primaryVertex(nullptr);
459 ATH_CHECK( defaultPVContainer.isValid() );
460 if (defaultPVContainer.cptr()->size()==0) {
461 ATH_MSG_WARNING("You have no primary vertices: " << defaultPVContainer.cptr()->size());
462 return StatusCode::RECOVERABLE;
463 }
464 else primaryVertex = (*defaultPVContainer.cptr())[0];
465
466 //----------------------------------------------------
467 // Record refitted primary vertices
468 //----------------------------------------------------
470 if(m_refitPV) {
472 ATH_CHECK( refPvContainer.record(std::make_unique<xAOD::VertexContainer>(), std::make_unique<xAOD::VertexAuxContainer>()) );
473 }
474
475 // Get TrackParticle container (standard + LRT)
477 ATH_CHECK( trackContainer.isValid() );
478
479 // Get all track containers when m_RelinkContainers is not empty
480 std::vector<const xAOD::TrackParticleContainer*> trackCols;
483 ATH_CHECK( handle.isValid() );
484 trackCols.push_back(handle.cptr());
485 }
486
487 // output V0 vertices
489 if(m_vertexV0ContainerKey.key()=="" && m_v0VtxOutputKey.key()!="") {
491 ATH_CHECK( V0OutputContainer.record(std::make_unique<xAOD::VertexContainer>(), std::make_unique<xAOD::VertexAuxContainer>()) );
492 }
493
494 // Get the input containers
495 // Note: If the event does not contain a JX candidate, it is skipped and the V0 container will not be constructed if not previously available in StoreGate
497 ATH_CHECK( jxContainer.isValid() );
498 if(jxContainer->size()==0) return StatusCode::SUCCESS;
499
500 // Select the displaced tracks
501 std::vector<const xAOD::TrackParticle*> tracksDisplaced;
502 if(m_v0VtxOutputKey.key()!="") {
503 for(const xAOD::TrackParticle* TP : *trackContainer.cptr()) {
504 // V0 track selection (https://gitlab.cern.ch/atlas/athena/-/blob/main/InnerDetector/InDetRecTools/InDetTrackSelectorTool/src/InDetConversionTrackSelectorTool.cxx)
505 if(m_v0TrkSelector->decision(*TP, primaryVertex)) {
506 uint8_t temp(0);
507 uint8_t nclus(0);
508 if(TP->summaryValue(temp, xAOD::numberOfPixelHits)) nclus += temp;
509 if(TP->summaryValue(temp, xAOD::numberOfSCTHits) ) nclus += temp;
510 if(!m_useTRT && nclus == 0) continue;
511
512 bool trk_cut = false;
513 if(nclus != 0) trk_cut = true;
514 if(nclus == 0 && TP->pt()>=m_ptTRT) trk_cut = true;
515 if(!trk_cut) continue;
516
517 // track is used if std::abs(d0/sig_d0) > d0_cut for PV
518 if(!d0Pass(TP,primaryVertex)) continue;
519
520 tracksDisplaced.push_back(TP);
521 }
522 }
523 }
524
525 SG::Accessor<std::string> mAcc_type("Type_V0Vtx");
526 SG::Accessor<int> mAcc_gfit("gamma_fit");
527 SG::Accessor<float> mAcc_gmass("gamma_mass");
528 SG::Accessor<float> mAcc_gchisq("gamma_chisq");
529 SG::Accessor<int> mAcc_gndof("gamma_ndof");
530
531 std::vector<std::pair<const xAOD::Vertex*,V0Enum> > selectedV0Candidates;
532
534 if(m_vertexV0ContainerKey.key() != "") {
536 ATH_CHECK( V0Container.isValid() );
537
538 for(const xAOD::Vertex* vtx : *V0Container.cptr()) {
539 std::string type_V0Vtx;
540 if(mAcc_type.isAvailable(*vtx)) type_V0Vtx = mAcc_type(*vtx);
541
542 V0Enum opt(UNKNOWN); double massV0(0);
543 if(type_V0Vtx == "Lambda") {
544 opt = LAMBDA;
545 massV0 = m_V0Tools->invariantMass(vtx, m_massesV0_ppi);
546 if(massV0<m_LambdaMassLower || massV0>m_LambdaMassUpper) continue;
547 }
548 else if(type_V0Vtx == "Lambdabar") {
549 opt = LAMBDABAR;
550 massV0 = m_V0Tools->invariantMass(vtx, m_massesV0_pip);
551 if(massV0<m_LambdaMassLower || massV0>m_LambdaMassUpper) continue;
552 }
553 else if(type_V0Vtx == "Ks") {
554 opt = KS;
555 massV0 = m_V0Tools->invariantMass(vtx, m_massesV0_pipi);
556 if(massV0<m_KsMassLower || massV0>m_KsMassUpper) continue;
557 }
558
559 if(opt==UNKNOWN) continue;
561 if((opt==LAMBDA || opt==LAMBDABAR) && m_V01Hypothesis == "Ks") continue;
562 if(opt==KS && m_V01Hypothesis == "Lambda") continue;
563 }
564
565 int gamma_fit = mAcc_gfit.isAvailable(*vtx) ? mAcc_gfit(*vtx) : 0;
566 double gamma_mass = mAcc_gmass.isAvailable(*vtx) ? mAcc_gmass(*vtx) : -1;
567 double gamma_chisq = mAcc_gchisq.isAvailable(*vtx) ? mAcc_gchisq(*vtx) : 999999;
568 double gamma_ndof = mAcc_gndof.isAvailable(*vtx) ? mAcc_gndof(*vtx) : 0;
569 double chisqPerNdof = (gamma_ndof > 0.) ? gamma_chisq/gamma_ndof : -999999.;
570 if(gamma_fit==1 && (gamma_mass < m_minMass_gamma) && (chisqPerNdof < m_chi2cut_gamma)) continue;
571
572 selectedV0Candidates.push_back(std::pair<const xAOD::Vertex*,V0Enum>{vtx,opt});
573 }
574 }
575 else {
576 // fit V0 vertices
577 fitV0Container(V0OutputContainer.ptr(), tracksDisplaced, trackCols);
578
579 for(const xAOD::Vertex* vtx : *V0OutputContainer.cptr()) {
580 std::string type_V0Vtx;
581 if(mAcc_type.isAvailable(*vtx)) type_V0Vtx = mAcc_type(*vtx);
582
583 V0Enum opt(UNKNOWN); double massV0(0);
584 if(type_V0Vtx == "Lambda") {
585 opt = LAMBDA;
586 massV0 = m_V0Tools->invariantMass(vtx, m_massesV0_ppi);
587 if(massV0<m_LambdaMassLower || massV0>m_LambdaMassUpper) continue;
588 }
589 else if(type_V0Vtx == "Lambdabar") {
590 opt = LAMBDABAR;
591 massV0 = m_V0Tools->invariantMass(vtx, m_massesV0_pip);
592 if(massV0<m_LambdaMassLower || massV0>m_LambdaMassUpper) continue;
593 }
594 else if(type_V0Vtx == "Ks") {
595 opt = KS;
596 massV0 = m_V0Tools->invariantMass(vtx, m_massesV0_pipi);
597 if(massV0<m_KsMassLower || massV0>m_KsMassUpper) continue;
598 }
599
600 if(opt==UNKNOWN) continue;
602 if((opt==LAMBDA || opt==LAMBDABAR) && m_V01Hypothesis == "Ks") continue;
603 if(opt==KS && m_V01Hypothesis == "Lambda") continue;
604 }
605
606 int gamma_fit = mAcc_gfit.isAvailable(*vtx) ? mAcc_gfit(*vtx) : 0;
607 double gamma_mass = mAcc_gmass.isAvailable(*vtx) ? mAcc_gmass(*vtx) : -1;
608 double gamma_chisq = mAcc_gchisq.isAvailable(*vtx) ? mAcc_gchisq(*vtx) : 999999;
609 double gamma_ndof = mAcc_gndof.isAvailable(*vtx) ? mAcc_gndof(*vtx) : 0;
610 if(gamma_fit==1 && gamma_mass<m_minMass_gamma && gamma_chisq/gamma_ndof<m_chi2cut_gamma) continue;
611
612 selectedV0Candidates.push_back(std::pair<const xAOD::Vertex*,V0Enum>{vtx,opt});
613 }
614 }
615
616 // sort and chop the V0 candidates
617 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(); } );
618 if(m_maxV0Candidates>0 && selectedV0Candidates.size()>m_maxV0Candidates) {
619 selectedV0Candidates.erase(selectedV0Candidates.begin()+m_maxV0Candidates, selectedV0Candidates.end());
620 }
621 if(selectedV0Candidates.size()==0) return StatusCode::SUCCESS;
622
623 std::vector<Trk::VxCascadeInfo*> cascadeinfoContainer;
624 ATH_CHECK( performSearch(cascadeinfoContainer, selectedV0Candidates, ctx) );
625
626 // sort and chop the main candidates
627 std::sort( cascadeinfoContainer.begin(), cascadeinfoContainer.end(), [](Trk::VxCascadeInfo* a, Trk::VxCascadeInfo* b) { return a->fitChi2()/a->nDoF() < b->fitChi2()/b->nDoF(); } );
628 if(m_maxMainVCandidates>0 && cascadeinfoContainer.size()>m_maxMainVCandidates) {
629 for(auto it=cascadeinfoContainer.begin()+m_maxMainVCandidates; it!=cascadeinfoContainer.end(); it++) delete *it;
630 cascadeinfoContainer.erase(cascadeinfoContainer.begin()+m_maxMainVCandidates, cascadeinfoContainer.end());
631 }
632
634 ATH_CHECK( evt.isValid() );
635 BPhysPVCascadeTools helper(&(*m_CascadeTools), evt.cptr());
636 helper.SetMinNTracksInPV(m_PV_minNTracks);
637
638 // Decorators for the main vertex: chi2, ndf, pt and pt error, plus the V0 vertex variables
639 SG::Decorator<VertexLinkVector> CascadeLinksDecor("CascadeVertexLinks");
640 SG::Decorator<float> chi2_decor("ChiSquared");
641 SG::Decorator<int> ndof_decor("nDoF");
642 SG::Decorator<float> Pt_decor("Pt");
643 SG::Decorator<float> PtErr_decor("PtErr");
644
645 SG::Decorator<float> lxy_SV1_decor("lxy_SV1");
646 SG::Decorator<float> lxyErr_SV1_decor("lxyErr_SV1");
647 SG::Decorator<float> a0xy_SV1_decor("a0xy_SV1");
648 SG::Decorator<float> a0xyErr_SV1_decor("a0xyErr_SV1");
649 SG::Decorator<float> a0z_SV1_decor("a0z_SV1");
650 SG::Decorator<float> a0zErr_SV1_decor("a0zErr_SV1");
651
652 SG::Decorator<float> lxy_SV2_decor("lxy_SV2");
653 SG::Decorator<float> lxyErr_SV2_decor("lxyErr_SV2");
654 SG::Decorator<float> a0xy_SV2_decor("a0xy_SV2");
655 SG::Decorator<float> a0xyErr_SV2_decor("a0xyErr_SV2");
656 SG::Decorator<float> a0z_SV2_decor("a0z_SV2");
657 SG::Decorator<float> a0zErr_SV2_decor("a0zErr_SV2");
658
659 SG::Decorator<float> lxy_SV3_decor("lxy_SV3");
660 SG::Decorator<float> lxyErr_SV3_decor("lxyErr_SV3");
661 SG::Decorator<float> a0xy_SV3_decor("a0xy_SV3");
662 SG::Decorator<float> a0xyErr_SV3_decor("a0xyErr_SV3");
663 SG::Decorator<float> a0z_SV3_decor("a0z_SV3");
664 SG::Decorator<float> a0zErr_SV3_decor("a0zErr_SV3");
665
666 SG::Decorator<float> chi2_V3_decor("ChiSquared_V3");
667 SG::Decorator<int> ndof_V3_decor("nDoF_V3");
668
669 for(auto cascade_info : cascadeinfoContainer) {
670 if(cascade_info==nullptr) {
671 ATH_MSG_ERROR("CascadeInfo is null");
672 continue;
673 }
674
675 const std::vector<xAOD::Vertex*> &cascadeVertices = cascade_info->vertices();
676 if(cascadeVertices.size() != topoN) ATH_MSG_ERROR("Incorrect number of vertices");
677 for(size_t i=0; i<topoN; i++) {
678 if(cascadeVertices[i]==nullptr) ATH_MSG_ERROR("Error null vertex");
679 }
680
681 cascade_info->setSVOwnership(false); // Prevent Container from deleting vertices
682 const auto mainVertex = cascadeVertices[topoN-1]; // this is the mother vertex
683 const std::vector< std::vector<TLorentzVector> > &moms = cascade_info->getParticleMoms();
684
685 // Identify the input JX
686 int ijx = m_JXSubVtx ? topoN-2 : topoN-1;
687 const xAOD::Vertex* jxVtx(nullptr);
688 if(m_jxDaug_num==4) jxVtx = FindVertex<4>(jxContainer.ptr(), cascadeVertices[ijx]);
689 else if(m_jxDaug_num==3) jxVtx = FindVertex<3>(jxContainer.ptr(), cascadeVertices[ijx]);
690 else jxVtx = FindVertex<2>(jxContainer.ptr(), cascadeVertices[ijx]);
691
692 xAOD::BPhysHypoHelper vtx(m_hypoName, mainVertex);
693
694 // Get refitted track momenta from all vertices, charged tracks only
695 BPhysPVCascadeTools::SetVectorInfo(vtx, cascade_info);
696 vtx.setPass(true);
697
698 //
699 // Decorate main vertex
700 //
701 // mass, mass error
702 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/VxCascadeInfo.h
703 BPHYS_CHECK( vtx.setMass(m_CascadeTools->invariantMass(moms[topoN-1])) );
704 BPHYS_CHECK( vtx.setMassErr(m_CascadeTools->invariantMassError(moms[topoN-1],cascade_info->getCovariance()[topoN-1])) );
705 // pt and pT error (the default pt of mainVertex is != the pt of the full cascade fit!)
706 Pt_decor(*mainVertex) = m_CascadeTools->pT(moms[topoN-1]);
707 PtErr_decor(*mainVertex) = m_CascadeTools->pTError(moms[topoN-1],cascade_info->getCovariance()[topoN-1]);
708 // chi2 and ndof (the default chi2 of mainVertex is != the chi2 of the full cascade fit!)
709 chi2_decor(*mainVertex) = cascade_info->fitChi2();
710 ndof_decor(*mainVertex) = cascade_info->nDoF();
711
712 // decorate the cascade vertices
713 lxy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],mainVertex);
714 lxyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->lxyError(moms[0],cascade_info->getCovariance()[0],cascadeVertices[0],mainVertex);
715 a0z_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0z(moms[0],cascadeVertices[0],mainVertex);
716 a0zErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0zError(moms[0],cascade_info->getCovariance()[0],cascadeVertices[0],mainVertex);
717 a0xy_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xy(moms[0],cascadeVertices[0],mainVertex);
718 a0xyErr_SV1_decor(*cascadeVertices[0]) = m_CascadeTools->a0xyError(moms[0],cascade_info->getCovariance()[0],cascadeVertices[0],mainVertex);
719
721 lxy_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[2]);
722 lxyErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->lxyError(moms[1],cascade_info->getCovariance()[1],cascadeVertices[1],cascadeVertices[2]);
723 a0z_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0z(moms[1],cascadeVertices[1],cascadeVertices[2]);
724 a0zErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0zError(moms[1],cascade_info->getCovariance()[1],cascadeVertices[1],cascadeVertices[2]);
725 a0xy_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0xy(moms[1],cascadeVertices[1],cascadeVertices[2]);
726 a0xyErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0xyError(moms[1],cascade_info->getCovariance()[1],cascadeVertices[1],cascadeVertices[2]);
727 }
728 else {
729 lxy_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->lxy(moms[1],cascadeVertices[1],mainVertex);
730 lxyErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->lxyError(moms[1],cascade_info->getCovariance()[1],cascadeVertices[1],mainVertex);
731 a0z_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0z(moms[1],cascadeVertices[1],mainVertex);
732 a0zErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0zError(moms[1],cascade_info->getCovariance()[1],cascadeVertices[1],mainVertex);
733 a0xy_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0xy(moms[1],cascadeVertices[1],mainVertex);
734 a0xyErr_SV2_decor(*cascadeVertices[1]) = m_CascadeTools->a0xyError(moms[1],cascade_info->getCovariance()[1],cascadeVertices[1],mainVertex);
735 }
736
737 if(m_JXSubVtx) {
738 lxy_SV3_decor(*cascadeVertices[2]) = m_CascadeTools->lxy(moms[2],cascadeVertices[2],mainVertex);
739 lxyErr_SV3_decor(*cascadeVertices[2]) = m_CascadeTools->lxyError(moms[2],cascade_info->getCovariance()[2],cascadeVertices[2],mainVertex);
740 a0z_SV3_decor(*cascadeVertices[2]) = m_CascadeTools->a0z(moms[2],cascadeVertices[2],mainVertex);
741 a0zErr_SV3_decor(*cascadeVertices[2]) = m_CascadeTools->a0zError(moms[2],cascade_info->getCovariance()[2],cascadeVertices[2],mainVertex);
742 a0xy_SV3_decor(*cascadeVertices[2]) = m_CascadeTools->a0xy(moms[2],cascadeVertices[2],mainVertex);
743 a0xyErr_SV3_decor(*cascadeVertices[2]) = m_CascadeTools->a0xyError(moms[2],cascade_info->getCovariance()[2],cascadeVertices[2],mainVertex);
744 }
745
746 chi2_V3_decor(*cascadeVertices[2]) = m_V0Tools->chisq(jxVtx);
747 ndof_V3_decor(*cascadeVertices[2]) = m_V0Tools->ndof(jxVtx);
748
749 if(m_cascadeFitWithPV==0) {
750 ATH_CHECK(helper.FillCandwithRefittedVertices(m_refitPV, defaultPVContainer.cptr(), m_refitPV ? refPvContainer.ptr() : 0, &(*m_pvRefitter), m_PV_max, m_DoVertexType, cascade_info, topoN-1, m_massMainV, vtx));
751 }
752
753 for(size_t i=0; i<topoN; i++) {
754 VtxWriteHandles[i].ptr()->push_back(cascadeVertices[i]);
755 }
756
757 // Set links to cascade vertices
758 VertexLinkVector cascadeVertexLinks;
759 VertexLink vertexLink1;
760 vertexLink1.setElement(cascadeVertices[0]);
761 vertexLink1.setStorableObject(*VtxWriteHandles[0].ptr());
762 if( vertexLink1.isValid() ) cascadeVertexLinks.push_back( vertexLink1 );
763 VertexLink vertexLink2;
764 vertexLink2.setElement(cascadeVertices[1]);
765 vertexLink2.setStorableObject(*VtxWriteHandles[1].ptr());
766 if( vertexLink2.isValid() ) cascadeVertexLinks.push_back( vertexLink2 );
767 if(topoN==4) {
768 VertexLink vertexLink3;
769 vertexLink3.setElement(cascadeVertices[2]);
770 vertexLink3.setStorableObject(*VtxWriteHandles[2].ptr());
771 if( vertexLink3.isValid() ) cascadeVertexLinks.push_back( vertexLink3 );
772 }
773 CascadeLinksDecor(*mainVertex) = std::move(cascadeVertexLinks);
774 } // loop over cascadeinfoContainer
775
776 // Deleting cascadeinfo since this won't be stored.
777 for (auto cascade_info : cascadeinfoContainer) {
778 if(cascade_info) delete cascade_info;
779 }
780
781 return StatusCode::SUCCESS;
782 }
783
784 bool JpsiXPlus2V0::d0Pass(const xAOD::TrackParticle* track, const xAOD::Vertex* PV) const {
785 bool pass = false;
786 const EventContext& ctx = Gaudi::Hive::currentContext();
787 std::unique_ptr<Trk::Perigee> per = m_trackToVertexTool->perigeeAtVertex(ctx, *track, PV->position());
788 if(!per) return pass;
789 double d0 = per->parameters()[Trk::d0];
790 double sig_d0 = sqrt((*per->covariance())(0,0));
791 if(std::abs(d0/sig_d0) > m_d0_cut) pass = true;
792 return pass;
793 }
794
795 Trk::VxCascadeInfo* JpsiXPlus2V0::fitMainVtx(const EventContext& ctx, const xAOD::Vertex* JXvtx, std::vector<double>& massesJX, const xAOD::Vertex* V01vtx, const V0Enum V01, const xAOD::Vertex* V02vtx, const V0Enum V02, const std::vector<const xAOD::TrackParticleContainer*>& trackCols, const xAOD::VertexContainer* defaultPVContainer, const xAOD::VertexContainer* pvContainer) const {
796 Trk::VxCascadeInfo* result(nullptr);
797
798 std::vector<const xAOD::TrackParticle*> tracksJX;
799 for(size_t i=0; i<JXvtx->nTrackParticles(); i++) tracksJX.push_back(JXvtx->trackParticle(i));
800 if (tracksJX.size() != massesJX.size()) {
801 ATH_MSG_ERROR("Problems with JX input: number of tracks or track mass inputs is not correct!");
802 return result;
803 }
804 std::vector<const xAOD::TrackParticle*> tracksV01;
805 for(size_t j=0; j<V01vtx->nTrackParticles(); j++) tracksV01.push_back(V01vtx->trackParticle(j));
806 std::vector<const xAOD::TrackParticle*> tracksV02;
807 for(size_t j=0; j<V02vtx->nTrackParticles(); j++) tracksV02.push_back(V02vtx->trackParticle(j));
808
809 std::vector<const xAOD::TrackParticle*> tracksJpsi{tracksJX[0], tracksJX[1]};
810 std::vector<const xAOD::TrackParticle*> tracksX;
811 if(m_jxDaug_num>=3) tracksX.push_back(tracksJX[2]);
812 if(m_jxDaug_num==4) tracksX.push_back(tracksJX[3]);
813
814 std::vector<double> massesV01;
815 if(V01==LAMBDA) massesV01 = m_massesV0_ppi;
816 else if(V01==LAMBDABAR) massesV01 = m_massesV0_pip;
817 else if(V01==KS) massesV01 = m_massesV0_pipi;
818 std::vector<double> massesV02;
819 if(V02==LAMBDA) massesV02 = m_massesV0_ppi;
820 else if(V02==LAMBDABAR) massesV02 = m_massesV0_pip;
821 else if(V02==KS) massesV02 = m_massesV0_pipi;
822
823 TLorentzVector p4_moth, p4_JX, p4_v01, p4_v02, tmp;
824 for(size_t it=0; it<JXvtx->nTrackParticles(); it++) {
825 tmp.SetPtEtaPhiM(JXvtx->trackParticle(it)->pt(), JXvtx->trackParticle(it)->eta(), JXvtx->trackParticle(it)->phi(), massesJX[it]);
826 p4_moth += tmp; p4_JX += tmp;
827 }
828 xAOD::BPhysHelper V01_helper(V01vtx);
829 for(int it=0; it<V01_helper.nRefTrks(); it++) {
830 p4_moth += V01_helper.refTrk(it,massesV01[it]);
831 p4_v01 += V01_helper.refTrk(it,massesV01[it]);
832 }
833 xAOD::BPhysHelper V02_helper(V02vtx);
834 for(int it=0; it<V02_helper.nRefTrks(); it++) {
835 p4_moth += V02_helper.refTrk(it,massesV02[it]);
836 p4_v02 += V02_helper.refTrk(it,massesV02[it]);
837 }
838 double main_mass = p4_moth.M();
840 if(m_jxDaug_num==2 && m_massJpsi>0) main_mass += - p4_JX.M() + m_massJpsi;
841 else if(m_jxDaug_num>=3 && m_massJX>0) main_mass += - p4_JX.M() + m_massJX;
842 if((V01==LAMBDA || V01==LAMBDABAR) && m_massLd>0) main_mass += - p4_v01.M() + m_massLd;
843 else if(V01==KS && m_massKs>0) main_mass += - p4_v01.M() + m_massKs;
844 if((V02==LAMBDA || V02==LAMBDABAR) && m_massLd>0) main_mass += - p4_v02.M() + m_massLd;
845 else if(V02==KS && m_massKs>0) main_mass += - p4_v02.M() + m_massKs;
846 }
847 if (main_mass < m_MassLower || main_mass > m_MassUpper) return result;
848
850 double JXV02_mass = (p4_JX+p4_v02).M();
852 if(m_jxDaug_num==2 && m_massJpsi>0) JXV02_mass += - p4_JX.M() + m_massJpsi;
853 else if(m_jxDaug_num>=3 && m_massJX>0) JXV02_mass += - p4_JX.M() + m_massJX;
854 if((V02==LAMBDA || V02==LAMBDABAR) && m_massLd>0) JXV02_mass += - p4_v02.M() + m_massLd;
855 else if(V02==KS && m_massKs>0) JXV02_mass += - p4_v02.M() + m_massKs;
856 }
857 if (JXV02_mass < m_JXV02MassLower || JXV02_mass > m_JXV02MassUpper) return result;
858 }
859
865 xAOD::BPhysHelper JX_helper(JXvtx);
866 const xAOD::Vertex* origPv_xAOD = (m_cascadeFitWithPV>=1 && m_cascadeFitWithPV<=4) ? JX_helper.origPv(pvtype) : nullptr;
867 const xAOD::Vertex* pv_xAOD = (m_cascadeFitWithPV>=1 && m_cascadeFitWithPV<=4) ? JX_helper.pv(pvtype) : nullptr;
868 std::unique_ptr<Trk::RecVertex> pv_AOD;
869 if(pv_xAOD) pv_AOD = std::make_unique<Trk::RecVertex>(pv_xAOD->position(),pv_xAOD->covariancePosition(),pv_xAOD->numberDoF(),pv_xAOD->chiSquared());
870
871 SG::Decorator<float> chi2_V1_decor("ChiSquared_V1");
872 SG::Decorator<int> ndof_V1_decor("nDoF_V1");
873 SG::Decorator<std::string> type_V1_decor("Type_V1");
874 SG::Decorator<float> chi2_V2_decor("ChiSquared_V2");
875 SG::Decorator<int> ndof_V2_decor("nDoF_V2");
876 SG::Decorator<std::string> type_V2_decor("Type_V2");
877
878 SG::Accessor<int> mAcc_gfit("gamma_fit");
879 SG::Accessor<float> mAcc_gmass("gamma_mass");
880 SG::Accessor<float> mAcc_gmasserr("gamma_massError");
881 SG::Accessor<float> mAcc_gchisq("gamma_chisq");
882 SG::Accessor<int> mAcc_gndof("gamma_ndof");
883 SG::Accessor<float> mAcc_gprob("gamma_probability");
884
885 SG::Decorator<int> mDec_gfit("gamma_fit");
886 SG::Decorator<float> mDec_gmass("gamma_mass");
887 SG::Decorator<float> mDec_gmasserr("gamma_massError");
888 SG::Decorator<float> mDec_gchisq("gamma_chisq");
889 SG::Decorator<int> mDec_gndof("gamma_ndof");
890 SG::Decorator<float> mDec_gprob("gamma_probability");
891 SG::Decorator< std::vector<float> > trk_pxDeco("TrackPx_V0nc");
892 SG::Decorator< std::vector<float> > trk_pyDeco("TrackPy_V0nc");
893 SG::Decorator< std::vector<float> > trk_pzDeco("TrackPz_V0nc");
894
895 std::vector<float> trk_px;
896 std::vector<float> trk_py;
897 std::vector<float> trk_pz;
898
899 // Apply the user's settings to the fitter
900 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
901 // Robustness: http://cdsweb.cern.ch/record/685551
902 int robustness = 0;
903 m_iVertexFitter->setRobustness(robustness, *state);
904 // Build up the topology
905 // Vertex list
906 std::vector<Trk::VertexID> vrtList;
907 // https://gitlab.cern.ch/atlas/athena/-/blob/main/Tracking/TrkVertexFitter/TrkVKalVrtFitter/TrkVKalVrtFitter/IVertexCascadeFitter.h
908 // V01 vertex
909 Trk::VertexID vID1;
910 if (m_constrV01) {
911 vID1 = m_iVertexFitter->startVertex(tracksV01,massesV01,*state,V01==KS?m_massKs:m_massLd);
912 } else {
913 vID1 = m_iVertexFitter->startVertex(tracksV01,massesV01,*state);
914 }
915 vrtList.push_back(vID1);
916 // V02 vertex
917 Trk::VertexID vID2;
918 if (m_constrV02) {
919 vID2 = m_iVertexFitter->nextVertex(tracksV02,massesV02,*state,V02==KS?m_massKs:m_massLd);
920 } else {
921 vID2 = m_iVertexFitter->nextVertex(tracksV02,massesV02,*state);
922 }
923 vrtList.push_back(vID2);
924 Trk::VertexID vID3;
925 if(m_JXSubVtx) {
926 if(m_JXV02SubVtx) { // for e.g. Lambda_b -> Jpsi + Lambda
927 // JX+V02 vertex
928 std::vector<Trk::VertexID> vrtList1{vID1};
929 std::vector<Trk::VertexID> vrtList2{vID2};
930 if (m_constrJXV02) {
931 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList2,*state,m_massJXV02);
932 } else {
933 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList2,*state);
934 }
935 vrtList1.push_back(vID3);
936 if (m_constrJX && m_jxDaug_num>2) {
937 std::vector<Trk::VertexID> cnstV; cnstV.clear();
938 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
939 ATH_MSG_WARNING("addMassConstraint for JX failed");
940 }
941 }
942 // Mother vertex including JX and two V0's
943 std::vector<const xAOD::TrackParticle*> tp; tp.clear();
944 std::vector<double> tp_masses; tp_masses.clear();
945 if(m_constrMainV) {
946 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList1,*state,m_massMainV);
947 } else {
948 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList1,*state);
949 }
950 }
951 else { // no JXV02SubVtx
952 // JX vertex
953 if (m_constrJX && m_jxDaug_num>2) {
954 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state,m_massJX);
955 } else {
956 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,*state);
957 }
958 vrtList.push_back(vID3);
959 // Mother vertex including JX and two V0's
960 std::vector<const xAOD::TrackParticle*> tp; tp.clear();
961 std::vector<double> tp_masses; tp_masses.clear();
962 if(m_constrMainV) {
963 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList,*state,m_massMainV);
964 } else {
965 m_iVertexFitter->nextVertex(tp,tp_masses,vrtList,*state);
966 }
967 }
968 }
969 else { // m_JXSubVtx=false
970 // Mother vertex including JX and two V0's
971 if(m_constrMainV) {
972 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state,m_massMainV);
973 } else {
974 vID3 = m_iVertexFitter->nextVertex(tracksJX,massesJX,vrtList,*state);
975 }
976 if (m_constrJX && m_jxDaug_num>2) {
977 std::vector<Trk::VertexID> cnstV; cnstV.clear();
978 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJX,cnstV,*state,m_massJX).isSuccess() ) {
979 ATH_MSG_WARNING("addMassConstraint for JX failed");
980 }
981 }
982 }
983
984 if (m_constrJpsi) {
985 std::vector<Trk::VertexID> cnstV; cnstV.clear();
986 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksJpsi,cnstV,*state,m_massJpsi).isSuccess() ) {
987 ATH_MSG_WARNING("addMassConstraint for Jpsi failed");
988 }
989 }
990 if (m_constrX && m_jxDaug_num==4 && m_massX>0) {
991 std::vector<Trk::VertexID> cnstV; cnstV.clear();
992 if ( !m_iVertexFitter->addMassConstraint(vID3,tracksX,cnstV,*state,m_massX).isSuccess() ) {
993 ATH_MSG_WARNING("addMassConstraint for X failed");
994 }
995 }
996 // Do the work
997 std::unique_ptr<Trk::VxCascadeInfo> fit_result = std::unique_ptr<Trk::VxCascadeInfo>( m_iVertexFitter->fitCascade(*state, pv_AOD.get(), pv_AOD.get() && m_firstDecayAtPV ? true : false) );
998
999 if (fit_result != nullptr) {
1000 for(auto v : fit_result->vertices()) {
1001 if(v->nTrackParticles()==0) {
1002 std::vector<ElementLink<xAOD::TrackParticleContainer> > nullLinkVector;
1003 v->setTrackParticleLinks(nullLinkVector);
1004 }
1005 }
1006 // reset links to original tracks
1007 BPhysPVCascadeTools::PrepareVertexLinks(fit_result.get(), trackCols);
1008
1009 // necessary to prevent memory leak
1010 fit_result->setSVOwnership(true);
1011
1012 // Chi2/DOF cut
1013 double chi2DOF = fit_result->fitChi2()/fit_result->nDoF();
1014 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
1015
1016 const std::vector<std::vector<TLorentzVector> > &moms = fit_result->getParticleMoms();
1017 const std::vector<xAOD::Vertex*> &cascadeVertices = fit_result->vertices();
1018 size_t iMoth = cascadeVertices.size()-1;
1019 double lxy_SV1 = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[iMoth]);
1020 double lxy_SV2 = (m_JXSubVtx && m_JXV02SubVtx) ? m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[2]) : m_CascadeTools->lxy(moms[1],cascadeVertices[1],cascadeVertices[iMoth]);
1021 if(chi2CutPassed && lxy_SV1>m_lxyV01_cut && lxy_SV2>m_lxyV02_cut) {
1022 chi2_V1_decor(*cascadeVertices[0]) = V01vtx->chiSquared();
1023 ndof_V1_decor(*cascadeVertices[0]) = V01vtx->numberDoF();
1024 if(V01==LAMBDA) type_V1_decor(*cascadeVertices[0]) = "Lambda";
1025 else if(V01==LAMBDABAR) type_V1_decor(*cascadeVertices[0]) = "Lambdabar";
1026 else if(V01==KS) type_V1_decor(*cascadeVertices[0]) = "Ks";
1027 mDec_gfit(*cascadeVertices[0]) = mAcc_gfit.isAvailable(*V01vtx) ? mAcc_gfit(*V01vtx) : 0;
1028 mDec_gmass(*cascadeVertices[0]) = mAcc_gmass.isAvailable(*V01vtx) ? mAcc_gmass(*V01vtx) : -1;
1029 mDec_gmasserr(*cascadeVertices[0]) = mAcc_gmasserr.isAvailable(*V01vtx) ? mAcc_gmasserr(*V01vtx) : -1;
1030 mDec_gchisq(*cascadeVertices[0]) = mAcc_gchisq.isAvailable(*V01vtx) ? mAcc_gchisq(*V01vtx) : 999999;
1031 mDec_gndof(*cascadeVertices[0]) = mAcc_gndof.isAvailable(*V01vtx) ? mAcc_gndof(*V01vtx) : 0;
1032 mDec_gprob(*cascadeVertices[0]) = mAcc_gprob.isAvailable(*V01vtx) ? mAcc_gprob(*V01vtx) : -1;
1033 trk_px.clear(); trk_py.clear(); trk_pz.clear();
1034 for(int it=0; it<V01_helper.nRefTrks(); it++) {
1035 trk_px.push_back( V01_helper.refTrk(it).Px() );
1036 trk_py.push_back( V01_helper.refTrk(it).Py() );
1037 trk_pz.push_back( V01_helper.refTrk(it).Pz() );
1038 }
1039 trk_pxDeco(*cascadeVertices[0]) = trk_px;
1040 trk_pyDeco(*cascadeVertices[0]) = trk_py;
1041 trk_pzDeco(*cascadeVertices[0]) = trk_pz;
1042
1043 chi2_V2_decor(*cascadeVertices[1]) = V02vtx->chiSquared();
1044 ndof_V2_decor(*cascadeVertices[1]) = V02vtx->numberDoF();
1045 if(V02==LAMBDA) type_V2_decor(*cascadeVertices[1]) = "Lambda";
1046 else if(V02==LAMBDABAR) type_V2_decor(*cascadeVertices[1]) = "Lambdabar";
1047 else if(V02==KS) type_V2_decor(*cascadeVertices[1]) = "Ks";
1048 mDec_gfit(*cascadeVertices[1]) = mAcc_gfit.isAvailable(*V02vtx) ? mAcc_gfit(*V02vtx) : 0;
1049 mDec_gmass(*cascadeVertices[1]) = mAcc_gmass.isAvailable(*V02vtx) ? mAcc_gmass(*V02vtx) : -1;
1050 mDec_gmasserr(*cascadeVertices[1]) = mAcc_gmasserr.isAvailable(*V02vtx) ? mAcc_gmasserr(*V02vtx) : -1;
1051 mDec_gchisq(*cascadeVertices[1]) = mAcc_gchisq.isAvailable(*V02vtx) ? mAcc_gchisq(*V02vtx) : 999999;
1052 mDec_gndof(*cascadeVertices[1]) = mAcc_gndof.isAvailable(*V02vtx) ? mAcc_gndof(*V02vtx) : 0;
1053 mDec_gprob(*cascadeVertices[1]) = mAcc_gprob.isAvailable(*V02vtx) ? mAcc_gprob(*V02vtx) : -1;
1054 trk_px.clear(); trk_py.clear(); trk_pz.clear();
1055 for(int it=0; it<V02_helper.nRefTrks(); it++) {
1056 trk_px.push_back( V02_helper.refTrk(it).Px() );
1057 trk_py.push_back( V02_helper.refTrk(it).Py() );
1058 trk_pz.push_back( V02_helper.refTrk(it).Pz() );
1059 }
1060 trk_pxDeco(*cascadeVertices[1]) = std::move(trk_px);
1061 trk_pyDeco(*cascadeVertices[1]) = std::move(trk_py);
1062 trk_pzDeco(*cascadeVertices[1]) = std::move(trk_pz);
1063
1064 result = fit_result.release();
1065 }
1066 }
1067
1068 if(pv_xAOD && result && result->getParticleMoms().size()>0) {
1069 size_t index = result->getParticleMoms().size() - 1;
1070 const std::vector<TLorentzVector> &mom = result->getParticleMoms()[index];
1071 const Amg::MatrixX &cov = result->getCovariance()[index];
1072 const xAOD::Vertex* mainVertex = result->vertices()[index];
1073 xAOD::BPhysHypoHelper vtx(m_hypoName, mainVertex);
1074 bool isInDefaultPVCont = false;
1075 for(const xAOD::Vertex* pvVtx : *defaultPVContainer) {
1076 if(pv_xAOD == pvVtx) { isInDefaultPVCont = true; break; }
1077 }
1078 if(isInDefaultPVCont) vtx.setPv( pv_xAOD, defaultPVContainer, pvtype );
1079 else vtx.setPv( pv_xAOD, pvContainer, pvtype );
1080 if(origPv_xAOD) vtx.setOrigPv( origPv_xAOD, defaultPVContainer, pvtype );
1081 vtx.setLxy ( m_CascadeTools->lxy (mom, vtx.vtx(), pv_xAOD), pvtype );
1082 vtx.setLxyErr ( m_CascadeTools->lxyError (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
1083 vtx.setA0 ( m_CascadeTools->a0 (mom, vtx.vtx(), pv_xAOD), pvtype );
1084 vtx.setA0Err ( m_CascadeTools->a0Error (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
1085 vtx.setA0xy ( m_CascadeTools->a0xy (mom, vtx.vtx(), pv_xAOD), pvtype );
1086 vtx.setA0xyErr( m_CascadeTools->a0xyError(mom, cov, vtx.vtx(), pv_xAOD), pvtype );
1087 vtx.setZ0 ( m_CascadeTools->a0z (mom, vtx.vtx(), pv_xAOD), pvtype );
1088 vtx.setZ0Err ( m_CascadeTools->a0zError (mom, cov, vtx.vtx(), pv_xAOD), pvtype );
1089 vtx.setRefitPVStatus( 0, pvtype );
1090 // Proper decay times
1091 vtx.setTau( m_CascadeTools->tau(mom, vtx.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
1092 vtx.setTauErr( m_CascadeTools->tauError(mom, cov, vtx.vtx(), pv_xAOD), pvtype, xAOD::BPhysHypoHelper::TAU_INV_MASS );
1093 vtx.setTau( m_CascadeTools->tau(mom, vtx.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
1094 vtx.setTauErr( m_CascadeTools->tauError(mom, cov, vtx.vtx(), pv_xAOD, m_massMainV), pvtype, xAOD::BPhysHypoHelper::TAU_CONST_MASS );
1095 }
1096
1097 return result;
1098 }
1099
1100 void JpsiXPlus2V0::fitV0Container(xAOD::VertexContainer* V0ContainerNew, const std::vector<const xAOD::TrackParticle*>& selectedTracks, const std::vector<const xAOD::TrackParticleContainer*>& trackCols) const {
1101 const EventContext& ctx = Gaudi::Hive::currentContext();
1102
1103 SG::Decorator<std::string> mDec_type("Type_V0Vtx");
1104 SG::Decorator<int> mDec_gfit("gamma_fit");
1105 SG::Decorator<float> mDec_gmass("gamma_mass");
1106 SG::Decorator<float> mDec_gmasserr("gamma_massError");
1107 SG::Decorator<float> mDec_gchisq("gamma_chisq");
1108 SG::Decorator<int> mDec_gndof("gamma_ndof");
1109 SG::Decorator<float> mDec_gprob("gamma_probability");
1110
1111 std::vector<const xAOD::TrackParticle*> posTracks;
1112 std::vector<const xAOD::TrackParticle*> negTracks;
1113 for(const xAOD::TrackParticle* TP : selectedTracks) {
1114 if(TP->charge()>0) posTracks.push_back(TP);
1115 else negTracks.push_back(TP);
1116 }
1117
1118 for(const xAOD::TrackParticle* TP1 : posTracks) {
1119 const Trk::Perigee& aPerigee1 = TP1->perigeeParameters();
1120 for(const xAOD::TrackParticle* TP2 : negTracks) {
1121 const Trk::Perigee& aPerigee2 = TP2->perigeeParameters();
1122 int sflag(0), errorcode(0);
1123 Amg::Vector3D startingPoint = m_vertexEstimator->getCirclesIntersectionPoint(&aPerigee1,&aPerigee2,sflag,errorcode);
1124 if (errorcode != 0) {startingPoint(0) = 0.0; startingPoint(1) = 0.0; startingPoint(2) = 0.0;}
1125
1126 if (errorcode == 0 || errorcode == 5 || errorcode == 6 || errorcode == 8) {
1127 Trk::PerigeeSurface perigeeSurface(startingPoint);
1128 const Trk::TrackParameters* extrapolatedPerigee1 = m_extrapolator->extrapolate(ctx,TP1->perigeeParameters(), perigeeSurface).release();
1129 const Trk::TrackParameters* extrapolatedPerigee2 = m_extrapolator->extrapolate(ctx,TP2->perigeeParameters(), perigeeSurface).release();
1130 std::vector<std::unique_ptr<const Trk::TrackParameters> > cleanup;
1131 if(!extrapolatedPerigee1) extrapolatedPerigee1 = &TP1->perigeeParameters();
1132 else cleanup.push_back(std::unique_ptr<const Trk::TrackParameters>(extrapolatedPerigee1));
1133 if(!extrapolatedPerigee2) extrapolatedPerigee2 = &TP2->perigeeParameters();
1134 else cleanup.push_back(std::unique_ptr<const Trk::TrackParameters>(extrapolatedPerigee2));
1135 if(extrapolatedPerigee1 && extrapolatedPerigee2) {
1136 bool pass = false;
1137 TLorentzVector v1; TLorentzVector v2;
1138 if(!pass) {
1139 v1.SetXYZM(extrapolatedPerigee1->momentum().x(),extrapolatedPerigee1->momentum().y(),extrapolatedPerigee1->momentum().z(),m_mass_proton);
1140 v2.SetXYZM(extrapolatedPerigee2->momentum().x(),extrapolatedPerigee2->momentum().y(),extrapolatedPerigee2->momentum().z(),m_mass_pion);
1141 if((v1+v2).M()>1030.0 && (v1+v2).M()<1200.0) pass = true;
1142 }
1143 if(!pass) {
1144 v1.SetXYZM(extrapolatedPerigee1->momentum().x(),extrapolatedPerigee1->momentum().y(),extrapolatedPerigee1->momentum().z(),m_mass_pion);
1145 v2.SetXYZM(extrapolatedPerigee2->momentum().x(),extrapolatedPerigee2->momentum().y(),extrapolatedPerigee2->momentum().z(),m_mass_proton);
1146 if((v1+v2).M()>1030.0 && (v1+v2).M()<1200.0) pass = true;
1147 }
1148 if(!pass) {
1149 v1.SetXYZM(extrapolatedPerigee1->momentum().x(),extrapolatedPerigee1->momentum().y(),extrapolatedPerigee1->momentum().z(),m_mass_pion);
1150 v2.SetXYZM(extrapolatedPerigee2->momentum().x(),extrapolatedPerigee2->momentum().y(),extrapolatedPerigee2->momentum().z(),m_mass_pion);
1151 if((v1+v2).M()>430.0 && (v1+v2).M()<565.0) pass = true;
1152 }
1153 if(pass) {
1154 std::vector<const xAOD::TrackParticle*> tracksV0;
1155 tracksV0.push_back(TP1); tracksV0.push_back(TP2);
1156 std::unique_ptr<xAOD::Vertex> V0vtx = m_iV0Fitter->fit(ctx, tracksV0, startingPoint);
1157 if(V0vtx && V0vtx->chiSquared()>=0) {
1158 double chi2DOF = V0vtx->chiSquared()/V0vtx->numberDoF();
1159 if(chi2DOF>m_chi2cut_V0) continue;
1160
1161 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);
1162 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);
1163 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);
1164 if(massSig_V0_Lambda1<=massSig_V0_Lambda2 && massSig_V0_Lambda1<=massSig_V0_Ks) {
1165 mDec_type(*V0vtx.get()) = "Lambda";
1166 }
1167 else if(massSig_V0_Lambda2<=massSig_V0_Lambda1 && massSig_V0_Lambda2<=massSig_V0_Ks) {
1168 mDec_type(*V0vtx.get()) = "Lambdabar";
1169 }
1170 else if(massSig_V0_Ks<=massSig_V0_Lambda1 && massSig_V0_Ks<=massSig_V0_Lambda2) {
1171 mDec_type(*V0vtx.get()) = "Ks";
1172 }
1173
1174 int gamma_fit = 0; int gamma_ndof = 0; double gamma_chisq = 999999.;
1175 double gamma_prob = -1., gamma_mass = -1., gamma_massErr = -1.;
1176 std::unique_ptr<xAOD::Vertex> gammaVtx = m_iGammaFitter->fit(ctx, tracksV0, m_V0Tools->vtx(V0vtx.get()));
1177 if (gammaVtx) {
1178 gamma_fit = 1;
1179 gamma_mass = m_V0Tools->invariantMass(gammaVtx.get(),m_mass_e,m_mass_e);
1180 gamma_massErr = m_V0Tools->invariantMassError(gammaVtx.get(),m_mass_e,m_mass_e);
1181 gamma_chisq = m_V0Tools->chisq(gammaVtx.get());
1182 gamma_ndof = m_V0Tools->ndof(gammaVtx.get());
1183 gamma_prob = m_V0Tools->vertexProbability(gammaVtx.get());
1184 }
1185 mDec_gfit(*V0vtx.get()) = gamma_fit;
1186 mDec_gmass(*V0vtx.get()) = gamma_mass;
1187 mDec_gmasserr(*V0vtx.get()) = gamma_massErr;
1188 mDec_gchisq(*V0vtx.get()) = gamma_chisq;
1189 mDec_gndof(*V0vtx.get()) = gamma_ndof;
1190 mDec_gprob(*V0vtx.get()) = gamma_prob;
1191
1192 xAOD::BPhysHelper V0_helper(V0vtx.get());
1193 V0_helper.setRefTrks(); // AOD only method
1194
1195 if(not trackCols.empty()){
1196 try {
1197 JpsiUpsilonCommon::RelinkVertexTracks(trackCols, V0vtx.get());
1198 } catch (std::runtime_error const& e) {
1199 ATH_MSG_ERROR(e.what());
1200 return;
1201 }
1202 }
1203
1204 V0ContainerNew->push_back(std::move(V0vtx));
1205 }
1206 }
1207 }
1208 }
1209 }
1210 }
1211 }
1212
1213 template<size_t NTracks>
1215 for (const xAOD::Vertex* v1 : *cont) {
1216 assert(v1->nTrackParticles() == NTracks);
1217 std::array<const xAOD::TrackParticle*, NTracks> a1;
1218 std::array<const xAOD::TrackParticle*, NTracks> a2;
1219 for(size_t i=0; i<NTracks; i++){
1220 a1[i] = v1->trackParticle(i);
1221 a2[i] = v->trackParticle(i);
1222 }
1223 std::sort(a1.begin(), a1.end());
1224 std::sort(a2.begin(), a2.end());
1225 if(a1 == a2) return v1;
1226 }
1227 return nullptr;
1228 }
1229}
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_DEBUG(x,...)
#define ATH_MSG_ERROR(x,...)
#define ATH_MSG_WARNING(x,...)
#define ATH_MSG_FATAL(x,...)
#define BPHYS_CHECK(EXP)
Useful CHECK macro.
: B-physics xAOD helpers.
ATLAS-specific HepMC functions.
bool passed(DecisionID id, const DecisionIDContainer &)
checks if required decision ID is in the set of IDs in the container
static Double_t a
value_type push_back(value_type pElem)
Add an element to the end of the collection.
static void PrepareVertexLinks(Trk::VxCascadeInfo *result, const xAOD::TrackParticleContainer *importedTrackCollection)
static void SetVectorInfo(xAOD::BPhysHelper &, const Trk::VxCascadeInfo *)
static void RelinkVertexTracks(std::span< const xAOD::TrackParticleContainer *const > trkcols, xAOD::Vertex *vtx)
SG::ReadHandleKey< xAOD::VertexContainer > m_pvContainerName
std::vector< std::string > m_vertexJXHypoNames
JpsiXPlus2V0(const std::string &type, const std::string &name, const IInterface *parent)
std::vector< double > m_massesV0_pip
PublicToolHandle< Trk::V0Tools > m_V0Tools
PublicToolHandle< Analysis::PrimaryVertexRefitter > m_pvRefitter
virtual StatusCode addBranches(const EventContext &ctx) const override
ToolHandle< Reco::ITrackToVertex > m_trackToVertexTool
ToolHandle< Trk::ITrackSelectorTool > m_v0TrkSelector
std::vector< double > m_massesV0_ppi
SG::ReadHandleKey< xAOD::VertexContainer > m_VxPrimaryCandidateName
ToolHandle< Trk::IExtrapolator > m_extrapolator
ToolHandle< Trk::TrkVKalVrtFitter > m_iVertexFitter
SG::ReadHandleKey< xAOD::EventInfo > m_eventInfo_key
ToolHandle< InDet::VertexPointEstimator > m_vertexEstimator
std::vector< double > m_massesV0_pipi
ToolHandle< Trk::IVertexFitter > m_iGammaFitter
bool d0Pass(const xAOD::TrackParticle *track, const xAOD::Vertex *PV) const
SG::WriteHandleKey< xAOD::VertexContainer > m_refPVContainerName
SG::WriteHandleKeyArray< xAOD::VertexContainer > m_cascadeOutputKeys
Trk::VxCascadeInfo * fitMainVtx(const EventContext &ctx, const xAOD::Vertex *JXvtx, std::vector< double > &massesJX, const xAOD::Vertex *V01vtx, const V0Enum V01, const xAOD::Vertex *V02vtx, const V0Enum V02, const std::vector< const xAOD::TrackParticleContainer * > &trackCols, const xAOD::VertexContainer *defaultPVContainer, const xAOD::VertexContainer *pvContainer) const
virtual StatusCode initialize() override
ToolHandle< Trk::TrkV0VertexFitter > m_iV0Fitter
const xAOD::Vertex * FindVertex(const xAOD::VertexContainer *cont, const xAOD::Vertex *v) const
PublicToolHandle< DerivationFramework::CascadeTools > m_CascadeTools
SG::WriteHandleKey< xAOD::VertexContainer > m_v0VtxOutputKey
void fitV0Container(xAOD::VertexContainer *V0ContainerNew, const std::vector< const xAOD::TrackParticle * > &selectedTracks, const std::vector< const xAOD::TrackParticleContainer * > &trackCols) const
SG::ReadHandleKey< xAOD::VertexContainer > m_vertexJXContainerKey
SG::ReadHandleKey< xAOD::TrackParticleContainer > m_TrkParticleCollection
SG::ReadHandleKeyArray< xAOD::TrackParticleContainer > m_RelinkContainers
StatusCode performSearch(std::vector< Trk::VxCascadeInfo * > &cascadeinfoContainer, const std::vector< std::pair< const xAOD::Vertex *, V0Enum > > &selectedV0Candidates, const EventContext &ctx) const
SG::ReadHandleKey< xAOD::VertexContainer > m_vertexV0ContainerKey
Helper class to provide type-safe access to aux data.
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
Property holding a SG store/key/clid from which a ReadHandle is made.
const_pointer_type ptr()
Dereference the pointer.
virtual bool isValid() override final
Can the handle be successfully dereferenced?
const_pointer_type cptr()
Dereference the pointer.
Property holding a SG store/key/clid from which a WriteHandle is made.
const_pointer_type cptr() const
Dereference the pointer.
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
pointer_type ptr()
Dereference the pointer.
const Amg::Vector3D & momentum() const
Access method for the momentum.
Class describing the Line to which the Perigee refers to.
float setZ0(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
longitudinal impact parameter
const xAOD::Vertex * origPv(const pv_type vertexType=BPhysHelper::PV_MIN_A0)
original PV
bool setLxyErr(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
its error
bool setLxy(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
Set the transverse decay distance and its error measured between the refitted primary vertex of type ...
float setZ0Err(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
longitudinal impact parameter error
const xAOD::Vertex * vtx() const
Getter method for the cached xAOD::Vertex.
float setA0(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
Set the 3D and transverse impact parameters and their error.
bool setPv(const xAOD::Vertex *pv, const xAOD::VertexContainer *vertexContainer, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
Set the refitted collision vertex of type pv_type.
const xAOD::Vertex * pv(const pv_type vertexType=BPhysHelper::PV_MIN_A0)
Get the refitted collision vertex of type pv_type.
TVector3 refTrk(const size_t index)
Returns i-th refitted track 3-momentum.
pv_type
: Enum type of the PV
bool setRefTrks(std::vector< float > px, std::vector< float > py, std::vector< float > pz)
Sets refitted track momenta.
float setA0Err(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
3D impact parameter error
bool setRefitPVStatus(int code, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
Set the exitCode of the refitter for vertex of type pv_type.
float setA0xyErr(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
transverse impact parameter error
int nRefTrks()
Returns number of stored refitted track momenta.
float setA0xy(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
transverse impact parameter
bool setOrigPv(const xAOD::Vertex *pv, const xAOD::VertexContainer *vertexContainer, const pv_type vertexType=BPhysHelper::PV_MIN_A0)
Set the original collision vertex of type pv_type.
bool setMass(const float val)
Set given invariant mass and its error.
bool setPass(bool passVal)
get the pass flag for this hypothesis
bool setTau(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0, const tau_type tauType=BPhysHypoHelper::TAU_CONST_MASS)
: Set the proper decay time and error.
bool setMassErr(const float val)
invariant mass error
bool setTauErr(const float val, const pv_type vertexType=BPhysHelper::PV_MIN_A0, const tau_type tauType=BPhysHypoHelper::TAU_CONST_MASS)
proper decay time error
virtual double phi() const override final
The azimuthal angle ( ) of the particle (has range to .).
virtual double pt() const override final
The transverse momentum ( ) of the particle.
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
size_t nTrackParticles() const
Get the number of tracks associated with this vertex.
const TrackParticle * trackParticle(size_t i) const
Get the pointer to a given track that was used in vertex reco.
float numberDoF() const
Returns the number of degrees of freedom of the vertex fit as float.
float chiSquared() const
Returns the of the vertex fit as float.
const Amg::Vector3D & position() const
Returns the 3-pos.
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Eigen::Matrix< double, 3, 1 > Vector3D
THE reconstruction tool.
ElementLink< xAOD::VertexContainer > VertexLink
std::vector< VertexLink > VertexLinkVector
static const int PSI2S
static const int MUON
static const int BCPLUS
static const int ELECTRON
static const int K0S
static const int LAMBDAB0
static const int PIPLUS
static const int JPSI
static const int LAMBDA0
static const int PROTON
SG::Decorator< T, ALLOC > Decorator
Helper class to provide type-safe access to aux data, specialized for JaggedVecElt.
Definition AuxElement.h:576
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
@ d0
Definition ParamDefs.h:63
ParametersBase< TrackParametersDim, Charged > TrackParameters
Definition index.py:1
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
TrackParticle_v1 TrackParticle
Reference the current persistent version:
VertexContainer_v1 VertexContainer
Definition of the current "Vertex container version".
Vertex_v1 Vertex
Define the latest version of the vertex class.
@ numberOfSCTHits
number of hits in SCT [unit8_t].
@ numberOfPixelHits
these are the pixel hits, including the b-layer [unit8_t].