ATLAS Offline Software
Loading...
Searching...
No Matches
PsiPlusPsiSingleVertex.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3*/
5// PsiPlusPsiSingleVertex.cxx, (c) ATLAS Detector software
7#include "PsiPlusPsiSingleVertex.h"
8#include "TrkVertexFitterInterfaces/IVertexFitter.h"
9#include "TrkVKalVrtFitter/TrkVKalVrtFitter.h"
10#include "TrkVertexAnalysisUtils/V0Tools.h"
11#include "BPhysPVCascadeTools.h"
12#include "BPhysPVTools.h"
13#include "xAODTracking/VertexAuxContainer.h"
14#include "xAODBPhys/BPhysHypoHelper.h"
15#include "GeneratorModules/GenData.h"
16#include "TruthUtils/HepMCHelpers.h"
17#include <algorithm>
18
19namespace DerivationFramework {
20 typedef ElementLink<xAOD::VertexContainer> VertexLink;
21 typedef std::vector<VertexLink> VertexLinkVector;
22 typedef ElementLink<xAOD::TrackParticleContainer> TrackParticleLink;
23 typedef std::vector<TrackParticleLink> TrackParticleLinkVector;
24
25 PsiPlusPsiSingleVertex::PsiPlusPsiSingleVertex(const std::string& type, const std::string& name, const IInterface* parent) : base_class(type,name,parent),
26 m_vertexPsi1ContainerKey(""),
27 m_vertexPsi2ContainerKey(""),
28 m_outputsKeys({"Psi1Vtx", "Psi2Vtx", "MainVtx"}),
29 m_VxPrimaryCandidateName("PrimaryVertices"),
30 m_trackContainerName("InDetTrackParticles"),
31 m_eventInfo_key("EventInfo"),
32 m_jpsi1MassLower(0.0),
33 m_jpsi1MassUpper(20000.0),
34 m_diTrack1MassLower(-1.0),
35 m_diTrack1MassUpper(-1.0),
36 m_psi1MassLower(0.0),
37 m_psi1MassUpper(25000.0),
38 m_jpsi2MassLower(0.0),
39 m_jpsi2MassUpper(20000.0),
40 m_diTrack2MassLower(-1.0),
41 m_diTrack2MassUpper(-1.0),
42 m_psi2MassLower(0.0),
43 m_psi2MassUpper(25000.0),
44 m_MassLower(0.0),
45 m_MassUpper(31000.0),
46 m_vtx1Daug_num(4),
47 m_vtx1Daug1MassHypo(-1),
48 m_vtx1Daug2MassHypo(-1),
49 m_vtx1Daug3MassHypo(-1),
50 m_vtx1Daug4MassHypo(-1),
51 m_vtx2Daug_num(4),
52 m_vtx2Daug1MassHypo(-1),
53 m_vtx2Daug2MassHypo(-1),
54 m_vtx2Daug3MassHypo(-1),
55 m_vtx2Daug4MassHypo(-1),
56 m_maxCandidates(0),
57 m_mass_jpsi1(-1),
58 m_mass_diTrk1(-1),
59 m_mass_psi1(-1),
60 m_mass_jpsi2(-1),
61 m_mass_diTrk2(-1),
62 m_mass_psi2(-1),
63 m_constrJpsi1(false),
64 m_constrDiTrk1(false),
65 m_constrPsi1(false),
66 m_constrJpsi2(false),
67 m_constrDiTrk2(false),
68 m_constrPsi2(false),
69 m_chi2cut(-1.0),
70 m_iVertexFitter("Trk::TrkVKalVrtFitter"),
71 m_pvRefitter("Analysis::PrimaryVertexRefitter"),
72 m_V0Tools("Trk::V0Tools")
73 {
74 declareProperty("Psi1Vertices", m_vertexPsi1ContainerKey);
75 declareProperty("Psi2Vertices", m_vertexPsi2ContainerKey);
76 declareProperty("Psi1VtxHypoNames", m_vertexPsi1HypoNames);
77 declareProperty("Psi2VtxHypoNames", m_vertexPsi2HypoNames);
78 declareProperty("VxPrimaryCandidateName", m_VxPrimaryCandidateName);
79 declareProperty("TrackContainerName", m_trackContainerName);
80 declareProperty("RefPVContainerName", m_refPVContainerName = "RefittedPrimaryVertices");
81 declareProperty("Jpsi1MassLowerCut", m_jpsi1MassLower);
82 declareProperty("Jpsi1MassUpperCut", m_jpsi1MassUpper);
83 declareProperty("DiTrack1MassLower", m_diTrack1MassLower); // only effective when m_vtx1Daug_num=4
84 declareProperty("DiTrack1MassUpper", m_diTrack1MassUpper); // only effective when m_vtx1Daug_num=4
85 declareProperty("Psi1MassLowerCut", m_psi1MassLower);
86 declareProperty("Psi1MassUpperCut", m_psi1MassUpper);
87 declareProperty("Jpsi2MassLowerCut", m_jpsi2MassLower);
88 declareProperty("Jpsi2MassUpperCut", m_jpsi2MassUpper);
89 declareProperty("DiTrack2MassLower", m_diTrack2MassLower); // only effective when m_vtx2Daug_num=4
90 declareProperty("DiTrack2MassUpper", m_diTrack2MassUpper); // only effective when m_vtx2Daug_num=4
91 declareProperty("Psi2MassLowerCut", m_psi2MassLower);
92 declareProperty("Psi2MassUpperCut", m_psi2MassUpper);
93 declareProperty("MassLowerCut", m_MassLower);
94 declareProperty("MassUpperCut", m_MassUpper);
95 declareProperty("HypothesisName", m_hypoName = "TQ");
96 declareProperty("NumberOfPsi1Daughters", m_vtx1Daug_num);
97 declareProperty("Vtx1Daug1MassHypo", m_vtx1Daug1MassHypo);
98 declareProperty("Vtx1Daug2MassHypo", m_vtx1Daug2MassHypo);
99 declareProperty("Vtx1Daug3MassHypo", m_vtx1Daug3MassHypo);
100 declareProperty("Vtx1Daug4MassHypo", m_vtx1Daug4MassHypo);
101 declareProperty("NumberOfPsi2Daughters", m_vtx2Daug_num);
102 declareProperty("Vtx2Daug1MassHypo", m_vtx2Daug1MassHypo);
103 declareProperty("Vtx2Daug2MassHypo", m_vtx2Daug2MassHypo);
104 declareProperty("Vtx2Daug3MassHypo", m_vtx2Daug3MassHypo);
105 declareProperty("Vtx2Daug4MassHypo", m_vtx2Daug4MassHypo);
106 declareProperty("MaxCandidates", m_maxCandidates);
107 declareProperty("Jpsi1Mass", m_mass_jpsi1);
108 declareProperty("DiTrack1Mass", m_mass_diTrk1);
109 declareProperty("Psi1Mass", m_mass_psi1);
110 declareProperty("Jpsi2Mass", m_mass_jpsi2);
111 declareProperty("DiTrack2Mass", m_mass_diTrk2);
112 declareProperty("Psi2Mass", m_mass_psi2);
113 declareProperty("ApplyJpsi1MassConstraint", m_constrJpsi1);
114 declareProperty("ApplyDiTrk1MassConstraint", m_constrDiTrk1); // only effective when m_vtx1Daug_num=4
115 declareProperty("ApplyPsi1MassConstraint", m_constrPsi1);
116 declareProperty("ApplyJpsi2MassConstraint", m_constrJpsi2);
117 declareProperty("ApplyDiTrk2MassConstraint", m_constrDiTrk2); // only effective when m_vtx2Daug_num=4
118 declareProperty("ApplyPsi2MassConstraint", m_constrPsi2);
119 declareProperty("Chi2Cut", m_chi2cut);
120 declareProperty("RefitPV", m_refitPV = true);
121 declareProperty("MaxnPV", m_PV_max = 1000);
122 declareProperty("MinNTracksInPV", m_PV_minNTracks = 0);
123 declareProperty("DoVertexType", m_DoVertexType = 7);
124 declareProperty("TrkVertexFitterTool", m_iVertexFitter);
125 declareProperty("PVRefitter", m_pvRefitter);
126 declareProperty("V0Tools", m_V0Tools);
127 declareProperty("OutputVertexCollections", m_outputsKeys);
128 }
129
130 StatusCode PsiPlusPsiSingleVertex::initialize() {
131 // retrieving vertex Fitter
132 ATH_CHECK( m_iVertexFitter.retrieve() );
133
134 // retrieve PV refitter
135 ATH_CHECK( m_pvRefitter.retrieve() );
136
137 // retrieving the V0 tools
138 ATH_CHECK( m_V0Tools.retrieve() );
139
140 ATH_CHECK( m_vertexPsi1ContainerKey.initialize() );
141 ATH_CHECK( m_vertexPsi2ContainerKey.initialize() );
142 ATH_CHECK( m_VxPrimaryCandidateName.initialize() );
143 ATH_CHECK( m_trackContainerName.initialize() );
144 ATH_CHECK( m_refPVContainerName.initialize() );
145 ATH_CHECK( m_outputsKeys.initialize() );
146 ATH_CHECK( m_eventInfo_key.initialize() );
147
148 auto gendata = std::make_shared<GenData>();
149
150 if(m_mass_jpsi1 < 0.) m_mass_jpsi1 = gendata->particleMass(MC::JPSI).value();
151 if(m_mass_jpsi2 < 0.) m_mass_jpsi2 = gendata->particleMass(MC::JPSI).value();
152 if(m_mass_psi1 < 0.) m_mass_psi1 = gendata->particleMass(MC::PSI2S).value();
153 if(m_mass_psi2 < 0.) m_mass_psi2 = gendata->particleMass(MC::PSI2S).value();
154
155 if(m_vtx1Daug1MassHypo < 0.) m_vtx1Daug1MassHypo = gendata->particleMass(MC::MUON).value();
156 if(m_vtx1Daug2MassHypo < 0.) m_vtx1Daug2MassHypo = gendata->particleMass(MC::MUON).value();
157 if(m_vtx1Daug_num>=3 && m_vtx1Daug3MassHypo < 0.) m_vtx1Daug3MassHypo = gendata->particleMass(MC::PIPLUS).value();
158 if(m_vtx1Daug_num==4 && m_vtx1Daug4MassHypo < 0.) m_vtx1Daug4MassHypo = gendata->particleMass(MC::PIPLUS).value();
159 if(m_vtx2Daug1MassHypo < 0.) m_vtx2Daug1MassHypo = gendata->particleMass(MC::MUON).value();
160 if(m_vtx2Daug2MassHypo < 0.) m_vtx2Daug2MassHypo = gendata->particleMass(MC::MUON).value();
161 if(m_vtx2Daug_num>=3 && m_vtx2Daug3MassHypo < 0.) m_vtx2Daug3MassHypo = gendata->particleMass(MC::PIPLUS).value();
162 if(m_vtx2Daug_num==4 && m_vtx2Daug4MassHypo < 0.) m_vtx2Daug4MassHypo = gendata->particleMass(MC::PIPLUS).value();
163
164 return StatusCode::SUCCESS;
165 }
166
167 StatusCode PsiPlusPsiSingleVertex::addBranches(const EventContext& ctx) const {
168 if (m_vtx1Daug_num < 2 || m_vtx1Daug_num > 4 || m_vtx2Daug_num < 2 || m_vtx2Daug_num > 4) {
169 ATH_MSG_FATAL("Incorrect number of Psi daughters");
170 return StatusCode::FAILURE;
171 }
172
173 constexpr int topoN = 3;
174 if(m_outputsKeys.size() != topoN) {
175 ATH_MSG_FATAL("Incorrect number of VtxContainers");
176 return StatusCode::FAILURE;
177 }
178 std::array<SG::WriteHandle<xAOD::VertexContainer>, topoN> VtxWriteHandles; int ikey(0);
179 for(const SG::WriteHandleKey<xAOD::VertexContainer>& key : m_outputsKeys) {
180 VtxWriteHandles[ikey] = SG::WriteHandle<xAOD::VertexContainer>(key, ctx);
181 ATH_CHECK( VtxWriteHandles[ikey].record(std::make_unique<xAOD::VertexContainer>(), std::make_unique<xAOD::VertexAuxContainer>()) );
182 ikey++;
183 }
184
185 //----------------------------------------------------
186 // retrieve primary vertices
187 //----------------------------------------------------
188 SG::ReadHandle<xAOD::VertexContainer> pvContainer(m_VxPrimaryCandidateName, ctx);
189 ATH_CHECK( pvContainer.isValid() );
190 if (pvContainer.cptr()->size()==0) {
191 ATH_MSG_WARNING("You have no primary vertices: " << pvContainer.cptr()->size());
192 return StatusCode::RECOVERABLE;
193 }
194
195 // Get TrackParticle container (for setting links to the original tracks)
196 SG::ReadHandle<xAOD::TrackParticleContainer> trackContainer(m_trackContainerName, ctx);
197 ATH_CHECK( trackContainer.isValid() );
198
199 std::vector<const xAOD::TrackParticle*> tracksJpsi1;
200 std::vector<const xAOD::TrackParticle*> tracksDiTrk1;
201 std::vector<const xAOD::TrackParticle*> tracksPsi1;
202 std::vector<const xAOD::TrackParticle*> tracksJpsi2;
203 std::vector<const xAOD::TrackParticle*> tracksDiTrk2;
204 std::vector<const xAOD::TrackParticle*> tracksPsi2;
205 std::vector<const xAOD::TrackParticle*> inputTracks;
206 std::vector<double> massesPsi1;
207 massesPsi1.push_back(m_vtx1Daug1MassHypo);
208 massesPsi1.push_back(m_vtx1Daug2MassHypo);
209 if(m_vtx1Daug_num>=3) massesPsi1.push_back(m_vtx1Daug3MassHypo);
210 if(m_vtx1Daug_num==4) massesPsi1.push_back(m_vtx1Daug4MassHypo);
211 std::vector<double> massesPsi2;
212 massesPsi2.push_back(m_vtx2Daug1MassHypo);
213 massesPsi2.push_back(m_vtx2Daug2MassHypo);
214 if(m_vtx2Daug_num>=3) massesPsi2.push_back(m_vtx2Daug3MassHypo);
215 if(m_vtx2Daug_num==4) massesPsi2.push_back(m_vtx2Daug4MassHypo);
216 std::vector<double> massesInputTracks;
217 massesInputTracks.reserve(massesPsi1.size() + massesPsi2.size());
218 for(auto mass : massesPsi1) massesInputTracks.push_back(mass);
219 for(auto mass : massesPsi2) massesInputTracks.push_back(mass);
220
221 // Get Psi1 container
222 SG::ReadHandle<xAOD::VertexContainer> psi1Container(m_vertexPsi1ContainerKey, ctx);
223 ATH_CHECK( psi1Container.isValid() );
224 SG::ReadHandle<xAOD::VertexContainer> psi2Container(m_vertexPsi2ContainerKey, ctx);
225 ATH_CHECK( psi2Container.isValid() );
226
227 // Select the psi1 candidates before calling fit
228 std::vector<const xAOD::Vertex*> selectedPsi1Candidates;
229 for(auto vxcItr=psi1Container->cbegin(); vxcItr!=psi1Container->cend(); ++vxcItr) {
230 // Check the passed flag first
231 const xAOD::Vertex* vtx = *vxcItr;
232 bool passed = false;
233 for(size_t i=0; i<m_vertexPsi1HypoNames.size(); i++) {
234 SG::AuxElement::Accessor<Char_t> flagAcc("passed_"+m_vertexPsi1HypoNames[i]);
235 if(flagAcc.isAvailable(*vtx) && flagAcc(*vtx)) {
236 passed = true;
237 }
238 }
239 if(m_vertexPsi1HypoNames.size() && !passed) continue;
240
241 // Check psi1 candidate invariant mass and skip if need be
242 if(m_vtx1Daug_num>2) {
243 double mass_psi1 = m_V0Tools->invariantMass(*vxcItr, massesPsi1);
244 if (mass_psi1 < m_psi1MassLower || mass_psi1 > m_psi1MassUpper) continue;
245 }
246
247 TLorentzVector p4_mu1, p4_mu2;
248 p4_mu1.SetPtEtaPhiM( (*vxcItr)->trackParticle(0)->pt(),
249 (*vxcItr)->trackParticle(0)->eta(),
250 (*vxcItr)->trackParticle(0)->phi(), m_vtx1Daug1MassHypo);
251 p4_mu2.SetPtEtaPhiM( (*vxcItr)->trackParticle(1)->pt(),
252 (*vxcItr)->trackParticle(1)->eta(),
253 (*vxcItr)->trackParticle(1)->phi(), m_vtx1Daug2MassHypo);
254 double mass_jpsi1 = (p4_mu1 + p4_mu2).M();
255 if (mass_jpsi1 < m_jpsi1MassLower || mass_jpsi1 > m_jpsi1MassUpper) continue;
256
257 if(m_vtx1Daug_num==4 && m_diTrack1MassLower>=0 && m_diTrack1MassUpper>m_diTrack1MassLower) {
258 TLorentzVector p4_trk1, p4_trk2;
259 p4_trk1.SetPtEtaPhiM( (*vxcItr)->trackParticle(2)->pt(),
260 (*vxcItr)->trackParticle(2)->eta(),
261 (*vxcItr)->trackParticle(2)->phi(), m_vtx1Daug3MassHypo);
262 p4_trk2.SetPtEtaPhiM( (*vxcItr)->trackParticle(3)->pt(),
263 (*vxcItr)->trackParticle(3)->eta(),
264 (*vxcItr)->trackParticle(3)->phi(), m_vtx1Daug4MassHypo);
265 double mass_diTrk1 = (p4_trk1 + p4_trk2).M();
266 if (mass_diTrk1 < m_diTrack1MassLower || mass_diTrk1 > m_diTrack1MassUpper) continue;
267 }
268
269 selectedPsi1Candidates.push_back(*vxcItr);
270 }
271
272 // Select the psi2 candidates before calling fit
273 std::vector<const xAOD::Vertex*> selectedPsi2Candidates;
274 for(auto vxcItr=psi2Container->cbegin(); vxcItr!=psi2Container->cend(); ++vxcItr) {
275 // Check the passed flag first
276 const xAOD::Vertex* vtx = *vxcItr;
277 bool passed = false;
278 for(size_t i=0; i<m_vertexPsi2HypoNames.size(); i++) {
279 SG::AuxElement::Accessor<Char_t> flagAcc("passed_"+m_vertexPsi2HypoNames[i]);
280 if(flagAcc.isAvailable(*vtx) && flagAcc(*vtx)) {
281 passed = true;
282 }
283 }
284 if(m_vertexPsi2HypoNames.size() && !passed) continue;
285
286 // Check psi2 candidate invariant mass and skip if need be
287 if(m_vtx2Daug_num>2) {
288 double mass_psi2 = m_V0Tools->invariantMass(*vxcItr,massesPsi2);
289 if(mass_psi2 < m_psi2MassLower || mass_psi2 > m_psi2MassUpper) continue;
290 }
291
292 TLorentzVector p4_mu1, p4_mu2;
293 p4_mu1.SetPtEtaPhiM( (*vxcItr)->trackParticle(0)->pt(),
294 (*vxcItr)->trackParticle(0)->eta(),
295 (*vxcItr)->trackParticle(0)->phi(), m_vtx2Daug1MassHypo);
296 p4_mu2.SetPtEtaPhiM( (*vxcItr)->trackParticle(1)->pt(),
297 (*vxcItr)->trackParticle(1)->eta(),
298 (*vxcItr)->trackParticle(1)->phi(), m_vtx2Daug2MassHypo);
299 double mass_jpsi2 = (p4_mu1 + p4_mu2).M();
300 if (mass_jpsi2 < m_jpsi2MassLower || mass_jpsi2 > m_jpsi2MassUpper) continue;
301
302 if(m_vtx2Daug_num==4 && m_diTrack2MassLower>=0 && m_diTrack2MassUpper>m_diTrack2MassLower) {
303 TLorentzVector p4_trk1, p4_trk2;
304 p4_trk1.SetPtEtaPhiM( (*vxcItr)->trackParticle(2)->pt(),
305 (*vxcItr)->trackParticle(2)->eta(),
306 (*vxcItr)->trackParticle(2)->phi(), m_vtx2Daug3MassHypo);
307 p4_trk2.SetPtEtaPhiM( (*vxcItr)->trackParticle(3)->pt(),
308 (*vxcItr)->trackParticle(3)->eta(),
309 (*vxcItr)->trackParticle(3)->phi(), m_vtx2Daug4MassHypo);
310 double mass_diTrk2 = (p4_trk1 + p4_trk2).M();
311 if (mass_diTrk2 < m_diTrack2MassLower || mass_diTrk2 > m_diTrack2MassUpper) continue;
312 }
313 selectedPsi2Candidates.push_back(*vxcItr);
314 }
315
316 std::vector<std::pair<const xAOD::Vertex*, const xAOD::Vertex*> > candidatePairs;
317 for(auto psi1Itr=selectedPsi1Candidates.cbegin(); psi1Itr!=selectedPsi1Candidates.cend(); ++psi1Itr) {
318 tracksPsi1.clear();
319 tracksPsi1.reserve((*psi1Itr)->nTrackParticles());
320 for(size_t i=0; i<(*psi1Itr)->nTrackParticles(); i++) tracksPsi1.push_back((*psi1Itr)->trackParticle(i));
321 for(auto psi2Itr=selectedPsi2Candidates.cbegin(); psi2Itr!=selectedPsi2Candidates.cend(); ++psi2Itr) {
322 bool skip = false;
323 for(size_t j=0; j<(*psi2Itr)->nTrackParticles(); j++) {
324 if(std::find(tracksPsi1.cbegin(), tracksPsi1.cend(), (*psi2Itr)->trackParticle(j)) != tracksPsi1.cend()) { skip = true; break; }
325 }
326 if(skip) continue;
327 if(m_vertexPsi1ContainerKey.key() == m_vertexPsi2ContainerKey.key()) {
328 for(size_t ic=0; ic<candidatePairs.size(); ic++) {
329 const xAOD::Vertex* psi1Vertex = candidatePairs[ic].first;
330 const xAOD::Vertex* psi2Vertex = candidatePairs[ic].second;
331 if((psi1Vertex == *psi1Itr && psi2Vertex == *psi2Itr) || (psi1Vertex == *psi2Itr && psi2Vertex == *psi1Itr)) { skip = true; break; }
332 }
333 }
334 if(skip) continue;
335 candidatePairs.push_back(std::pair<const xAOD::Vertex*, const xAOD::Vertex*>(*psi1Itr,*psi2Itr));
336 }
337 }
338
339 std::sort( candidatePairs.begin(), candidatePairs.end(), [](std::pair<const xAOD::Vertex*, const xAOD::Vertex*> a, std::pair<const xAOD::Vertex*, const xAOD::Vertex*> b) { return a.first->chiSquared()/a.first->numberDoF()+a.second->chiSquared()/a.second->numberDoF() < b.first->chiSquared()/b.first->numberDoF()+b.second->chiSquared()/b.second->numberDoF(); } );
340 if(m_maxCandidates>0 && candidatePairs.size()>m_maxCandidates) {
341 candidatePairs.erase(candidatePairs.begin()+m_maxCandidates, candidatePairs.end());
342 }
343
344 for(size_t ic=0; ic<candidatePairs.size(); ic++) {
345 const xAOD::Vertex* psi1Vertex = candidatePairs[ic].first;
346 const xAOD::Vertex* psi2Vertex = candidatePairs[ic].second;
347
348 tracksPsi1.clear();
349 tracksPsi1.reserve(psi1Vertex->nTrackParticles());
350 for(size_t it=0; it<psi1Vertex->nTrackParticles(); it++) tracksPsi1.push_back(psi1Vertex->trackParticle(it));
351 if (tracksPsi1.size() != massesPsi1.size()) {
352 ATH_MSG_ERROR("Problems with Psi1 input: number of tracks or track mass inputs is not correct!");
353 }
354 tracksPsi2.clear();
355 tracksPsi2.reserve(psi2Vertex->nTrackParticles());
356 for(size_t it=0; it<psi2Vertex->nTrackParticles(); it++) tracksPsi2.push_back(psi2Vertex->trackParticle(it));
357 if (tracksPsi2.size() != massesPsi2.size()) {
358 ATH_MSG_ERROR("Problems with Psi2 input: number of tracks or track mass inputs is not correct!");
359 }
360
361 tracksJpsi1.clear();
362 tracksJpsi1.push_back(psi1Vertex->trackParticle(0));
363 tracksJpsi1.push_back(psi1Vertex->trackParticle(1));
364 tracksDiTrk1.clear();
365 if(m_vtx1Daug_num==4) {
366 tracksDiTrk1.push_back(psi1Vertex->trackParticle(2));
367 tracksDiTrk1.push_back(psi1Vertex->trackParticle(3));
368 }
369 tracksJpsi2.clear();
370 tracksJpsi2.push_back(psi2Vertex->trackParticle(0));
371 tracksJpsi2.push_back(psi2Vertex->trackParticle(1));
372 tracksDiTrk2.clear();
373 if(m_vtx2Daug_num==4) {
374 tracksDiTrk2.push_back(psi2Vertex->trackParticle(2));
375 tracksDiTrk2.push_back(psi2Vertex->trackParticle(3));
376 }
377
378 TLorentzVector p4_moth;
379 TLorentzVector tmp;
380 for(size_t it=0; it<psi1Vertex->nTrackParticles(); it++) {
381 tmp.SetPtEtaPhiM(psi1Vertex->trackParticle(it)->pt(),psi1Vertex->trackParticle(it)->eta(),psi1Vertex->trackParticle(it)->phi(),massesPsi1[it]);
382 p4_moth += tmp;
383 }
384 for(size_t it=0; it<psi2Vertex->nTrackParticles(); it++) {
385 tmp.SetPtEtaPhiM(psi2Vertex->trackParticle(it)->pt(),psi2Vertex->trackParticle(it)->eta(),psi2Vertex->trackParticle(it)->phi(),massesPsi2[it]);
386 p4_moth += tmp;
387 }
388 if (p4_moth.M() < m_MassLower || p4_moth.M() > m_MassUpper) continue;
389
390 inputTracks.clear();
391 inputTracks.push_back(psi1Vertex->trackParticle(0));
392 inputTracks.push_back(psi1Vertex->trackParticle(1));
393 if(m_vtx1Daug_num>=3) inputTracks.push_back(psi1Vertex->trackParticle(2));
394 if(m_vtx1Daug_num==4) inputTracks.push_back(psi1Vertex->trackParticle(3));
395 inputTracks.push_back(psi2Vertex->trackParticle(0));
396 inputTracks.push_back(psi2Vertex->trackParticle(1));
397 if(m_vtx2Daug_num>=3) inputTracks.push_back(psi2Vertex->trackParticle(2));
398 if(m_vtx2Daug_num==4) inputTracks.push_back(psi2Vertex->trackParticle(3));
399
400 // start the fit
401 std::unique_ptr<Trk::IVKalState> state = m_iVertexFitter->makeState(ctx);
402 m_iVertexFitter->setMassInputParticles(massesInputTracks, *state);
403 if (m_constrJpsi1) {
404 m_iVertexFitter->setMassForConstraint(m_mass_jpsi1, std::array<int,2>{1,2}, *state);
405 }
406 if (m_constrPsi1 && m_vtx1Daug_num>=3) {
407 m_iVertexFitter->setMassForConstraint(m_mass_psi1, m_vtx1Daug_num==4 ? std::vector<int>{1,2,3,4} : std::vector<int>{1,2,3}, *state);
408 }
409 if (m_constrDiTrk1 && m_vtx1Daug_num==4) {
410 m_iVertexFitter->setMassForConstraint(m_mass_diTrk1, std::array<int,2>{3,4}, *state);
411 }
412 if (m_constrJpsi2) {
413 m_iVertexFitter->setMassForConstraint(m_mass_jpsi2, std::array<int,2>{m_vtx1Daug_num+1,m_vtx1Daug_num+2}, *state);
414 }
415 if (m_constrPsi2 && m_vtx2Daug_num>=3) {
416 m_iVertexFitter->setMassForConstraint(m_mass_psi2, m_vtx1Daug_num==4 ? std::vector<int>{m_vtx1Daug_num+1,m_vtx1Daug_num+2,m_vtx1Daug_num+3,m_vtx1Daug_num+4} : std::vector<int>{m_vtx1Daug_num+1,m_vtx1Daug_num+2,m_vtx1Daug_num+3}, *state);
417 }
418 if (m_constrDiTrk2 && m_vtx2Daug_num==4) {
419 m_iVertexFitter->setMassForConstraint(m_mass_diTrk2, std::array<int,2>{m_vtx1Daug_num+3,m_vtx1Daug_num+4}, *state);
420 }
421
422 // Starting point
423 Amg::Vector3D startingPoint((psi1Vertex->x()+psi2Vertex->x())/2,(psi1Vertex->y()+psi2Vertex->y())/2,(psi1Vertex->z()+psi2Vertex->z())/2);
424
425 // do the fit
426 std::unique_ptr<xAOD::Vertex> theResult( m_iVertexFitter->fit(inputTracks, startingPoint, *state) );
427
428 if(theResult != nullptr){
429 // Chi2/DOF cut
430 double chi2DOF = theResult->chiSquared()/theResult->numberDoF();
431 bool chi2CutPassed = (m_chi2cut <= 0.0 || chi2DOF < m_chi2cut);
432 if(chi2CutPassed) {
433 TrackParticleLinkVector tpLinkVector;
434 for(size_t i=0; i<theResult->trackParticleLinks().size(); i++) {
435 TrackParticleLink mylink = theResult->trackParticleLinks()[i];
436 mylink.setStorableObject(*trackContainer.ptr(), true);
437 tpLinkVector.push_back( mylink );
438 }
439 theResult->clearTracks();
440 theResult->setTrackParticleLinks( tpLinkVector );
441
442 xAOD::Vertex* psi1Vertex_ = new xAOD::Vertex(*psi1Vertex);
443 xAOD::BPhysHelper psi1_helper(psi1Vertex_);
444 psi1_helper.setRefTrks();
445 TrackParticleLinkVector tpLinkVector_psi1;
446 for(size_t i=0; i<psi1Vertex_->trackParticleLinks().size(); i++) {
447 TrackParticleLink mylink = psi1Vertex_->trackParticleLinks()[i];
448 mylink.setStorableObject(*trackContainer.ptr(), true);
449 tpLinkVector_psi1.push_back( mylink );
450 }
451 psi1Vertex_->clearTracks();
452 psi1Vertex_->setTrackParticleLinks( tpLinkVector_psi1 );
453
454 xAOD::Vertex* psi2Vertex_ = new xAOD::Vertex(*psi2Vertex);
455 xAOD::BPhysHelper psi2_helper(psi2Vertex_);
456 psi2_helper.setRefTrks();
457 TrackParticleLinkVector tpLinkVector_psi2;
458 for(size_t i=0; i<psi2Vertex_->trackParticleLinks().size(); i++) {
459 TrackParticleLink mylink = psi2Vertex_->trackParticleLinks()[i];
460 mylink.setStorableObject(*trackContainer.ptr(), true);
461 tpLinkVector_psi2.push_back( mylink );
462 }
463 psi2Vertex_->clearTracks();
464 psi2Vertex_->setTrackParticleLinks( tpLinkVector_psi2 );
465
466 VtxWriteHandles[0].ptr()->push_back(psi1Vertex_);
467 VtxWriteHandles[1].ptr()->push_back(psi2Vertex_);
468
469 // Set links to preceding vertices
470 VertexLinkVector precedingVertexLinks;
471 VertexLink vertexLink1;
472 vertexLink1.setElement(psi1Vertex_);
473 vertexLink1.setStorableObject(*VtxWriteHandles[0].ptr());
474 if( vertexLink1.isValid() ) precedingVertexLinks.push_back( vertexLink1 );
475 VertexLink vertexLink2;
476 vertexLink2.setElement(psi2Vertex_);
477 vertexLink2.setStorableObject(*VtxWriteHandles[1].ptr());
478 if( vertexLink2.isValid() ) precedingVertexLinks.push_back( vertexLink2 );
479
480 SG::AuxElement::Decorator<VertexLinkVector> PrecedingLinksDecor("PrecedingVertexLinks");
481 PrecedingLinksDecor(*theResult.get()) = precedingVertexLinks;
482
483 xAOD::BPhysHypoHelper vtx(m_hypoName, theResult.get());
484 vtx.setRefTrks();
485 vtx.setPass(true);
486
487 VtxWriteHandles[2].ptr()->push_back(theResult.release()); //ptr is released preventing deletion
488 }
489 }
490 } //Iterate over candidatePairs vertices
491
492 SG::ReadHandle<xAOD::EventInfo> evt(m_eventInfo_key, ctx);
493 ATH_CHECK( evt.isValid() );
494 BPhysPVTools helper(&(*m_V0Tools), evt.cptr());
495 helper.SetMinNTracksInPV(m_PV_minNTracks);
496
497 if(m_refitPV) {
498 SG::WriteHandle<xAOD::VertexContainer> refPvContainer(m_refPVContainerName, ctx);
499 ATH_CHECK( refPvContainer.record(std::make_unique<xAOD::VertexContainer>(), std::make_unique<xAOD::VertexAuxContainer>()) );
500
501 if(VtxWriteHandles[2]->size()>0) {
502 for(int i=0; i<topoN; i++) {
503 StatusCode SC = helper.FillCandwithRefittedVertices(VtxWriteHandles[i].ptr(), pvContainer.cptr(), refPvContainer.ptr(), &(*m_pvRefitter), m_PV_max, m_DoVertexType);
504 if(SC.isFailure()){
505 ATH_MSG_FATAL("FillCandwithRefittedVertices failed - check the vertices you passed");
506 return SC;
507 }
508 }
509 }
510 }
511 else {
512 if(VtxWriteHandles[2]->size()>0) {
513 for(int i=0; i<topoN; i++) {
514 StatusCode SC = helper.FillCandExistingVertices(VtxWriteHandles[i].ptr(), pvContainer.cptr(), m_DoVertexType);
515 if(SC.isFailure()){
516 ATH_MSG_FATAL("FillCandExistingVertices failed - check the vertices you passed");
517 return SC;
518 }
519 }
520 }
521 }
522
523 return StatusCode::SUCCESS;
524 }
525}
ATH_CHECK
#define ATH_CHECK
Evaluate an expression and check for errors.
Definition AthCheckMacros.h:40
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition AthMsgStreamMacros.h:33
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition AthMsgStreamMacros.h:34
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition AthMsgStreamMacros.h:32
BPhysHypoHelper.h
: B-physics xAOD helpers.
HepMCHelpers.h
ATLAS-specific HepMC functions.
passed
bool passed(DecisionID id, const DecisionIDContainer &)
checks if required decision ID is in the set of IDs in the container
Definition TrigCompositeUtilsRoot.cxx:118
a
static Double_t a
Definition LArPhysWaveHECTool.cxx:38
size
size_t size() const
Number of registered mappings.
DerivationFramework::PsiPlusPsiSingleVertex::m_refPVContainerName
SG::WriteHandleKey< xAOD::VertexContainer > m_refPVContainerName
Definition PsiPlusPsiSingleVertex.h:86
DerivationFramework::PsiPlusPsiSingleVertex::addBranches
virtual StatusCode addBranches(const EventContext &ctx) const override
Definition PsiPlusPsiSingleVertex.cxx:167
DerivationFramework::PsiPlusPsiSingleVertex::m_trackContainerName
SG::ReadHandleKey< xAOD::TrackParticleContainer > m_trackContainerName
Definition PsiPlusPsiSingleVertex.h:38
DerivationFramework::PsiPlusPsiSingleVertex::m_eventInfo_key
SG::ReadHandleKey< xAOD::EventInfo > m_eventInfo_key
Definition PsiPlusPsiSingleVertex.h:39
DerivationFramework::PsiPlusPsiSingleVertex::m_VxPrimaryCandidateName
SG::ReadHandleKey< xAOD::VertexContainer > m_VxPrimaryCandidateName
Name of primary vertex container.
Definition PsiPlusPsiSingleVertex.h:37
DerivationFramework::PsiPlusPsiSingleVertex::m_pvRefitter
PublicToolHandle< Analysis::PrimaryVertexRefitter > m_pvRefitter
Definition PsiPlusPsiSingleVertex.h:82
DerivationFramework::PsiPlusPsiSingleVertex::m_outputsKeys
SG::WriteHandleKeyArray< xAOD::VertexContainer > m_outputsKeys
Definition PsiPlusPsiSingleVertex.h:36
DerivationFramework::PsiPlusPsiSingleVertex::m_vertexPsi2ContainerKey
SG::ReadHandleKey< xAOD::VertexContainer > m_vertexPsi2ContainerKey
Definition PsiPlusPsiSingleVertex.h:33
DerivationFramework::PsiPlusPsiSingleVertex::PsiPlusPsiSingleVertex
PsiPlusPsiSingleVertex(const std::string &t, const std::string &n, const IInterface *p)
Definition PsiPlusPsiSingleVertex.cxx:25
DerivationFramework::PsiPlusPsiSingleVertex::m_vertexPsi1ContainerKey
SG::ReadHandleKey< xAOD::VertexContainer > m_vertexPsi1ContainerKey
Definition PsiPlusPsiSingleVertex.h:32
DerivationFramework::PsiPlusPsiSingleVertex::m_iVertexFitter
ToolHandle< Trk::TrkVKalVrtFitter > m_iVertexFitter
Definition PsiPlusPsiSingleVertex.h:81
SG::ReadHandle::ptr
const_pointer_type ptr()
Dereference the pointer.
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
SG::ReadHandle::cptr
const_pointer_type cptr()
Dereference the pointer.
SG::WriteHandleKey
Property holding a SG store/key/clid from which a WriteHandle is made.
Definition StoreGate/StoreGate/WriteHandleKey.h:40
SG::WriteHandle::record
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
SG::WriteHandle::ptr
pointer_type ptr()
Dereference the pointer.
xAOD::BPhysHelper::setRefTrks
bool setRefTrks(std::vector< float > px, std::vector< float > py, std::vector< float > pz)
Sets refitted track momenta.
Definition BPhysHelper.cxx:286
xAOD::BPhysHypoHelper::setPass
bool setPass(bool passVal)
get the pass flag for this hypothesis
Definition BPhysHypoHelper.cxx:367
xAOD::TrackParticle_v1::phi
virtual double phi() const override final
The azimuthal angle ( ) of the particle (has range to .).
xAOD::TrackParticle_v1::pt
virtual double pt() const override final
The transverse momentum ( ) of the particle.
Definition TrackParticle_v1.cxx:75
xAOD::TrackParticle_v1::eta
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
Definition TrackParticle_v1.cxx:79
xAOD::Vertex_v1::z
float z() const
Returns the z position.
xAOD::Vertex_v1::setTrackParticleLinks
void setTrackParticleLinks(const TrackParticleLinks_t &trackParticles)
Set all track particle links at once.
xAOD::Vertex_v1::nTrackParticles
size_t nTrackParticles() const
Get the number of tracks associated with this vertex.
Definition Vertex_v1.cxx:292
xAOD::Vertex_v1::clearTracks
void clearTracks()
Remove all tracks from the vertex.
Definition Vertex_v1.cxx:353
xAOD::Vertex_v1::trackParticleLinks
const TrackParticleLinks_t & trackParticleLinks() const
Get all the particles associated with the vertex.
xAOD::Vertex_v1::trackParticle
const TrackParticle * trackParticle(size_t i) const
Get the pointer to a given track that was used in vertex reco.
Definition Vertex_v1.cxx:271
xAOD::Vertex_v1::y
float y() const
Returns the y position.
xAOD::Vertex_v1::x
float x() const
Returns the x position.
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47
DerivationFramework
THE reconstruction tool.
Definition CascadeTools.h:16
DerivationFramework::TrackParticleLink
ElementLink< xAOD::TrackParticleContainer > TrackParticleLink
Definition PsiPlusPsiSingleVertex.cxx:22
DerivationFramework::VertexLink
ElementLink< xAOD::VertexContainer > VertexLink
Definition JpsiPlusDpstCascade.cxx:23
DerivationFramework::VertexLinkVector
std::vector< VertexLink > VertexLinkVector
Definition JpsiPlusDpstCascade.cxx:24
DerivationFramework::TrackParticleLinkVector
std::vector< TrackParticleLink > TrackParticleLinkVector
Definition PsiPlusPsiSingleVertex.cxx:23
MC::PSI2S
static const int PSI2S
Definition HepMCHelpers.h:130
MC::MUON
static const int MUON
Definition HepMCHelpers.h:78
MC::PIPLUS
static const int PIPLUS
Definition HepMCHelpers.h:110
MC::JPSI
static const int JPSI
Definition HepMCHelpers.h:122
std::sort
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
Definition DVL_algorithms.h:554
xAOD::Vertex
Vertex_v1 Vertex
Define the latest version of the vertex class.
Definition Event/xAOD/xAODTracking/xAODTracking/Vertex.h:16
type
Generated on for ATLAS Offline Software by doxygen 1.16.1