ATLAS Offline Software
Loading...
Searching...
No Matches
MuPlusDsCascade.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3*/
5// MuPlusDsCascade.cxx, (c) ATLAS Detector software
7#include "MuPlusDsCascade.h"
8#include "DerivationFrameworkBPhys/CascadeTools.h"
9#include "BPhysPVCascadeTools.h"
10#include "xAODBPhys/BPhysHypoHelper.h"
11#include "TrkVertexFitterInterfaces/IVertexFitter.h"
12#include "TrkVKalVrtFitter/TrkVKalVrtFitter.h"
13#include "TrkVertexAnalysisUtils/V0Tools.h"
14#include "TrkToolInterfaces/ITrackSelectorTool.h"
15
16#include <algorithm>
17#include "HepPDT/ParticleDataTable.hh"
18
19#include "xAODMuon/MuonContainer.h"
20#include "xAODTracking/VertexContainer.h"
21#include "xAODTracking/VertexAuxContainer.h"
22#include "xAODEventInfo/EventInfo.h"
23#include <xAODCore/AuxContainerBase.h>
24
25#include "AthContainers/ConstDataVector.h"
26#include "TrkVKalVrtFitter/VxCascadeInfo.h"
27#include "InDetTrackSelectionTool/InDetTrackSelectionTool.h"
28#include "TruthUtils/HepMCHelpers.h"
29
30namespace DerivationFramework {
31 typedef ElementLink<xAOD::VertexContainer> VertexLink;
32 typedef std::vector<VertexLink> VertexLinkVector;
33 typedef std::vector<const xAOD::TrackParticle*> TrackBag;
34
35 typedef ElementLink<xAOD::MuonContainer> MuonsLink;
36 typedef std::vector<MuonsLink> MuonsLinkVector;
37
38 StatusCode MuPlusDsCascade::initialize() {
39
40 // retrieving vertex Fitter
41 ATH_CHECK( m_iVertexFitter.retrieve());
42
43 // retrieving the V0 tools //from JpsiFinder
44 ATH_CHECK( m_V0Tools.retrieve());
45
46 // Get the track selector tool from ToolSvc
47 if ( m_trkSelector.retrieve().isFailure() ) {
48 ATH_MSG_FATAL("Failed to retrieve tool " << m_trkSelector);
49 return StatusCode::FAILURE;
50 } else {
51 ATH_MSG_INFO("Retrieved tool " << m_trkSelector);
52 }
53
54 //======================== inDetTrack selection tool ==================
55 m_trackSelectionTools = std::make_unique<InDet::InDetTrackSelectionTool>("TrackSelector");
56 ANA_CHECK(m_trackSelectionTools->setProperty("CutLevel", "LoosePrimary"));
57 ANA_CHECK(m_trackSelectionTools->initialize() );
58 //=====================================================================
59
60 // retrieving the Cascade tools
61 ATH_CHECK( m_CascadeTools.retrieve());
62
63 ATH_CHECK( m_eventInfo_key.initialize() );
64
65 ATH_CHECK( m_partPropSvc.retrieve() );
66 auto pdt = m_partPropSvc->PDT();
67
68 // Ds+/- : K K π
69 // D+/- : K π π
70 // Lambda_c+/- : K π p
71
72 // retrieve particle masses
73 if(m_vtx0MassHypo < 0.)
74 m_vtx0MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::BCPLUS);
75 if(m_vtx1MassHypo < 0.) {
76 if(std::abs(m_Dx_pid) == 411) m_vtx1MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::DPLUS);
77 if(std::abs(m_Dx_pid) == 431) m_vtx1MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::DSPLUS);
78 if(std::abs(m_Dx_pid) == 4122) m_vtx1MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::LAMBDACPLUS);
79 }
80
81 if(m_vtx0Daug1MassHypo < 0.) m_vtx0Daug1MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::MUON);
82 if(m_vtx1Daug1MassHypo < 0.) {
83 if(std::abs(m_Dx_pid) == 431) m_vtx1Daug1MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::KPLUS); //Ds+
84 else m_vtx1Daug1MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::PIPLUS); //D+, Lambda_c+
85 }
86 if(m_vtx1Daug2MassHypo < 0.) m_vtx1Daug2MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::KPLUS);
87 if(m_vtx1Daug3MassHypo < 0.) {
88 if(std::abs(m_Dx_pid) == 4122) m_vtx1Daug3MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::PROTON); //Lambda_c+
89 else m_vtx1Daug3MassHypo = BPhysPVCascadeTools::getParticleMass(pdt, MC::PIPLUS); //Ds+, D+
90 }
91 return StatusCode::SUCCESS;
92 }
93
94
95 StatusCode MuPlusDsCascade::addBranches(const EventContext& ctx) const
96 {
97 std::vector<Trk::VxCascadeInfo*> cascadeinfoContainer;
98 constexpr int topoN = 2;
99 std::array<xAOD::VertexContainer*, topoN> Vtxwritehandles;
100 std::array<xAOD::VertexAuxContainer*, topoN> Vtxwritehandlesaux;
101 if(m_cascadeOutputsKeys.size() !=topoN) { ATH_MSG_FATAL("Incorrect number of VtxContainers"); return StatusCode::FAILURE; }
102
103 for(int i =0; i<topoN;i++){
104 Vtxwritehandles[i] = new xAOD::VertexContainer();
105 Vtxwritehandlesaux[i] = new xAOD::VertexAuxContainer();
106 Vtxwritehandles[i]->setStore(Vtxwritehandlesaux[i]);
107 ATH_CHECK(evtStore()->record(Vtxwritehandles[i] , m_cascadeOutputsKeys[i] )); // FIXME Use Handles
108 ATH_CHECK(evtStore()->record(Vtxwritehandlesaux[i], m_cascadeOutputsKeys[i] + "Aux.")); // FIXME Use Handles
109 }
110
111 //----------------------------------------------------
112 // retrieve primary vertices
113 //----------------------------------------------------
114 const xAOD::Vertex * primaryVertex{};
115 const xAOD::VertexContainer *pvContainer{};
116 ATH_CHECK(evtStore()->retrieve(pvContainer, m_VxPrimaryCandidateName)); // FIXME Use Handles
117 ATH_MSG_DEBUG("Found " << m_VxPrimaryCandidateName << " in StoreGate!");
118
119 if (pvContainer->size()==0){
120 ATH_MSG_WARNING("You have no primary vertices: " << pvContainer->size());
121 return StatusCode::RECOVERABLE;
122 } else {
123 primaryVertex = (*pvContainer)[0];
124 }
125
126 //----------------------------------------------------
127 // Try to retrieve refitted primary vertices
128 //----------------------------------------------------
129 xAOD::VertexContainer* refPvContainer{};
130 xAOD::VertexAuxContainer* refPvAuxContainer{};
131 if (m_refitPV) {
132 if (evtStore()->contains<xAOD::VertexContainer>(m_refPVContainerName)) {
133 // refitted PV container exists. Get it from the store gate
134 ATH_CHECK(evtStore()->retrieve(refPvContainer , m_refPVContainerName )); // FIXME Use Handles
135 ATH_CHECK(evtStore()->retrieve(refPvAuxContainer, m_refPVContainerName + "Aux.")); // FIXME Use Handles
136 } else {
137 // refitted PV container does not exist. Create a new one.
138 refPvContainer = new xAOD::VertexContainer;
139 refPvAuxContainer = new xAOD::VertexAuxContainer;
140 refPvContainer->setStore(refPvAuxContainer);
141 ATH_CHECK(evtStore()->record(refPvContainer , m_refPVContainerName)); // FIXME Use Handles
142 ATH_CHECK(evtStore()->record(refPvAuxContainer, m_refPVContainerName+"Aux.")); // FIXME Use Handles
143 }
144 }
145
146 ATH_CHECK(performSearch(&cascadeinfoContainer, ctx));
147 SG::ReadHandle<xAOD::EventInfo> evt(m_eventInfo_key, ctx);
148 if(!evt.isValid()) {
149 ATH_MSG_ERROR("Cannot Retrieve " << m_eventInfo_key.key() );
150 return StatusCode::FAILURE;
151 }
152
153 BPhysPVCascadeTools helper(&(*m_CascadeTools), evt.cptr());
154 helper.SetMinNTracksInPV(m_PV_minNTracks);
155
156 // Decorators for the main vertex: chi2, ndf, pt and pt error, plus the V0 vertex variables
157 SG::AuxElement::Decorator<VertexLinkVector> CascadeLinksDecor("CascadeVertexLinks");
158 SG::AuxElement::Decorator<MuonsLinkVector> MuonsLinksDecor("MuonsLinks");
159 SG::AuxElement::Decorator<VertexLinkVector> DxLinksDecor("DxVertexLinks");
160 SG::AuxElement::Decorator<float> chi2_decor("ChiSquared");
161 SG::AuxElement::Decorator<float> ndof_decor("NumberDoF");
162 SG::AuxElement::Decorator<float> Pt_decor("Pt");
163 SG::AuxElement::Decorator<float> PtErr_decor("PtErr");
164 SG::AuxElement::Decorator<float> Mass_svdecor("Dx_mass");
165 SG::AuxElement::Decorator<float> MassErr_svdecor("Dx_massErr");
166 SG::AuxElement::Decorator<float> Pt_svdecor("Dx_Pt");
167 SG::AuxElement::Decorator<float> PtErr_svdecor("Dx_PtErr");
168 SG::AuxElement::Decorator<float> Lxy_svdecor("Dx_Lxy");
169 SG::AuxElement::Decorator<float> LxyErr_svdecor("Dx_LxyErr");
170 SG::AuxElement::Decorator<float> Tau_svdecor("Dx_Tau");
171 SG::AuxElement::Decorator<float> TauErr_svdecor("Dx_TauErr");
172 SG::AuxElement::Decorator<float> Dx_chi2_svdecor("Dx_chi2");
173 SG::AuxElement::Decorator<float> Dx_ndof_svdecor("Dx_ndof");
174
175 // K X1 X2
176 // Ds+/- : K K π
177 // D+/- : K π π
178 // Lambda_c+/- : K π p
179
180 SG::AuxElement::Decorator<float> massKX1_svdecor("KX_mass");
181 SG::AuxElement::Decorator<float> massKX1X2_svdecor("KXpi_mass");
182 SG::AuxElement::Decorator<float> RapidityKX1X2_svdecor("KXpi_Rapidity");
183
184 SG::AuxElement::Decorator<float> MuMass_decor("Mu_mass");
185 SG::AuxElement::Decorator<float> MuPt_decor("Mu_pt");
186 SG::AuxElement::Decorator<float> MuEta_decor("Mu_eta");
187 SG::AuxElement::Decorator<float> MuChi2_decor("Mu_chi2");
188 SG::AuxElement::Decorator<float> MunDoF_decor("Mu_nDoF");
189 //muon contribution to chi2 of the cascade fit
190 SG::AuxElement::Decorator<float> MuChi2B_decor("Mu_chi2_B");
191 SG::AuxElement::Decorator<float> MunDoFB_decor("Mu_nDoF_B");
192
193 SG::AuxElement::Decorator<float> ChargeK_decor("K_charge");
194 SG::AuxElement::Decorator<float> ChargeX1_decor("X_charge"); // K (Ds+) or pi (D+, Lambda_c+)
195 SG::AuxElement::Decorator<float> ChargeX2_decor("Pi_charge"); // pi (Ds+, D+) or p (Lambda_c+)
196 SG::AuxElement::Decorator<float> ChargeMu_decor("Mu_charge");
197
198 ATH_MSG_DEBUG("cascadeinfoContainer size " << cascadeinfoContainer.size());
199
200 // Get Muons container
201 const xAOD::MuonContainer* muonContainer{};
202 ATH_CHECK( evtStore()->retrieve(muonContainer, m_muonCollectionKey) ); // FIXME Use Handles
203 ATH_MSG_DEBUG("Muon container size "<<muonContainer->size());
204
205 // Get D_(s)+/Lambda_c+ container and identify the input D_(s)+/Lambda_c+
206 const xAOD::VertexContainer *dxContainer{};
207 ATH_CHECK(evtStore()->retrieve(dxContainer , m_vertexDxContainerKey )); // FIXME Use Handles
208
209 for (Trk::VxCascadeInfo* x : cascadeinfoContainer) {
210 if(!x){
211 ATH_MSG_ERROR("cascadeinfoContainer is null");
212 continue;
213 }
214
215 // the cascade fitter returns:
216 // std::vector<xAOD::Vertex*>, each xAOD::Vertex contains the refitted track parameters (perigee at the vertex position)
217 // vertices[iv] the links to the original TPs and a covariance of size 3+5*NTRK; the chi2 of the total fit
218 // is split between the cascade vertices as per track contribution
219 // std::vector< std::vector<TLorentzVector> >, each std::vector<TLorentzVector> contains the refitted momenta (TLorentzVector)
220 // momenta[iv][...] of all tracks in the corresponding vertex, including any pseudotracks (from cascade vertices)
221 // originating in this vertex; the masses are as assigned in the cascade fit
222 // std::vector<Amg::MatrixX>, the corresponding covariance matrices in momentum space
223 // covariance[iv]
224 // int nDoF, double Chi2
225 //
226 // the invariant mass, pt, lifetime etc. errors should be calculated using the covariance matrices in momentum space as these
227 // take into account the full track-track and track-vertex correlations
228 //
229 // in the case of Jpsi+V0: vertices[0] is the V0 vertex, vertices[1] is the B/Lambda_b(bar) vertex, containing the 2 Jpsi tracks.
230 // The covariance terms between the two vertices are not stored. In momentum space momenta[0] contains the 2 V0 tracks,
231 // their momenta add up to the momentum of the 3rd track in momenta[1], the first two being the Jpsi tracks
232
233 const std::vector<xAOD::Vertex*> &cascadeVertices = x->vertices();
234
235 if(cascadeVertices.size()!=topoN) ATH_MSG_ERROR("Incorrect number of vertices");
236 if(cascadeVertices[0] == nullptr || cascadeVertices[1] == nullptr) ATH_MSG_ERROR("Error null vertex");
237 // Keep vertices (bear in mind that they come in reverse order!)
238 for(int i =0;i<topoN;i++) Vtxwritehandles[i]->push_back(cascadeVertices[i]);
239
240 x->setSVOwnership(false); // Prevent Container from deleting vertices
241 const auto mainVertex = cascadeVertices[1]; // this is the B_c+/- vertex
242 //and cascadeVertices[0] is Dx vertex
243 const std::vector< std::vector<TLorentzVector> > &moms = x->getParticleMoms();
244
245
246 // Set links to cascade vertices
247 std::vector<const xAOD::Vertex*> verticestoLink;
248 verticestoLink.push_back(cascadeVertices[0]);
249 if(Vtxwritehandles[1] == nullptr) ATH_MSG_ERROR("Vtxwritehandles[1] is null");
250 if(!BPhysPVCascadeTools::LinkVertices(CascadeLinksDecor, verticestoLink, Vtxwritehandles[0], cascadeVertices[1]))
251 ATH_MSG_ERROR("Error decorating with cascade vertices");
252
253 //----------------------------------- mu
254 // Identify muons
255 typedef std::vector<const xAOD::Muon*> MuonBag;
256 MuonBag selectedMuons; selectedMuons.clear();
257
258 for(const xAOD::Muon * mu : *muonContainer){
259 const xAOD::TrackParticle* muonTrk = mu->trackParticle(xAOD::Muon::InnerDetectorTrackParticle );
260 if (muonTrk == cascadeVertices[1]->trackParticle(0)) selectedMuons.push_back(mu); //there is always only one muon
261 }
262 ATH_MSG_DEBUG("selectedMuon size "<<selectedMuons.size()); //always only one muon
263
264 //-----------------------------------
265 // Create link for Muon
266 // create tmp vector of preceding vertex links
267 MuonsLinkVector preMuLinks;
268
269 // loop over input precedingMuons
270 auto muItr = selectedMuons.begin();
271 for(; muItr != selectedMuons.end(); muItr++) {
272 // create element link
273 MuonsLink muLink;
274 muLink.setElement(*muItr);
275 muLink.setStorableObject(*muonContainer);
276
277 // sanity check : is the link valid?
278 if( !muLink.isValid() ) continue;
279
280 // link is OK, store it in the tmp vector
281 preMuLinks.push_back( muLink );
282
283 } // end of loop over preceding vertices
284
285 // all OK: store preceding vertex links in the aux store
286 MuonsLinksDecor(*cascadeVertices[1]) = preMuLinks;
287 //-----------------------------------
288 //----------------------------------- mu
289
290 //----------------------------------- Dx / Lambda_c
291 // Identify the input D_(s)+ or Lambda_c+
292 const xAOD::Vertex* dxVertex = BPhysPVCascadeTools::FindVertex<3>(dxContainer, cascadeVertices[0]);
293 ATH_MSG_DEBUG("1 pt D_(s)+/Lambda_c+ tracks " << cascadeVertices[1]->trackParticle(0)->pt()<<", "
294 << cascadeVertices[0]->trackParticle(0)->pt()<< ", "<< cascadeVertices[0]->trackParticle(1)->pt()<< ", "
295 << cascadeVertices[0]->trackParticle(2)->pt());
296 if (dxVertex) ATH_MSG_DEBUG("2 pt D_(s)+/Lambda_c+ tracks " << dxVertex->trackParticle(0)->pt() << ", "
297 << dxVertex->trackParticle(1)->pt() << ", " << dxVertex->trackParticle(2)->pt());
298 // Set links to input vertices
299 std::vector<const xAOD::Vertex*> dxVerticestoLink;
300 if (dxVertex) dxVerticestoLink.push_back(dxVertex);
301 else ATH_MSG_WARNING("Could not find linking D_(s)+/Lambda_c+");
302 if(!BPhysPVCascadeTools::LinkVertices(DxLinksDecor, dxVerticestoLink, dxContainer, cascadeVertices[1]))
303 ATH_MSG_ERROR("Error decorating with D_(s)+/Lambda_c+ vertices");
304
305 // tag for D- and Lambda_c-
306 bool tagDp(true);
307 if (dxVertex) {
308 if((std::abs(m_Dx_pid)==411 || std::abs(m_Dx_pid)==4122) && (dxVertex->trackParticle(2)->charge()==-1)) tagDp = false;
309 }
310 //----------------------------------- Dx / Lambda_c
311
312 double mass_b = m_vtx0MassHypo;
313 double mass_d = m_vtx1MassHypo;
314 std::vector<double> massesMu;
315 massesMu.push_back(m_vtx0Daug1MassHypo); // µ
316 std::vector<double> massesDx;
317 if(tagDp){
318 massesDx.push_back(m_vtx1Daug1MassHypo); // K (Ds+/-) or π (D+, Lambda_c+)
319 massesDx.push_back(m_vtx1Daug2MassHypo); // K
320 }else{ // Change the order for D- or Lambda_c-
321 massesDx.push_back(m_vtx1Daug2MassHypo); // K
322 massesDx.push_back(m_vtx1Daug1MassHypo); // π
323 }
324 massesDx.push_back(m_vtx1Daug3MassHypo); // π (Ds+/-, D+-) or p (Lambda_c+)
325 std::vector<double> Masses;
326 Masses.push_back(m_vtx0Daug1MassHypo); // µ
327 Masses.push_back(m_vtx1MassHypo); // Dx / Lambda_c
328
329 // loop over candidates -- Don't apply PV_minNTracks requirement here
330 // because it may result in exclusion of the high-pt PV.
331 // get good PVs
332
333
334 xAOD::BPhysHypoHelper vtx(m_hypoName, mainVertex);
335
336 // Get refitted track momenta from all vertices, charged tracks only
337 BPhysPVCascadeTools::SetVectorInfo(vtx, x);
338
339 // Decorate main vertex
340 //
341 // 1.a) mass, mass error
342 BPHYS_CHECK( vtx.setMass(m_CascadeTools->invariantMass(moms[1])) );
343 BPHYS_CHECK( vtx.setMassErr(m_CascadeTools->invariantMassError(moms[1],x->getCovariance()[1])) );
344 // 1.b) pt and pT error (the default pt of mainVertex is != the pt of the full cascade fit!)
345 Pt_decor(*mainVertex) = m_CascadeTools->pT(moms[1]);
346 PtErr_decor(*mainVertex) = m_CascadeTools->pTError(moms[1],x->getCovariance()[1]);
347 // 1.c) chi2 and ndof (the default chi2 of mainVertex is != the chi2 of the full cascade fit!)
348 chi2_decor(*mainVertex) = x->fitChi2();
349 ndof_decor(*mainVertex) = x->nDoF();
350 Dx_chi2_svdecor(*mainVertex) = cascadeVertices[0]->chiSquared();
351 Dx_ndof_svdecor(*mainVertex) = cascadeVertices[0]->numberDoF();
352
353 //--------------------- mu
354 float muM = 0., mupt = 0., mueta = 0.;
355 if (selectedMuons.size()==1) {
356 TLorentzVector p4_mu1;
357 p4_mu1.SetPtEtaPhiM(selectedMuons.at(0)->pt(),
358 selectedMuons.at(0)->eta(),
359 selectedMuons.at(0)->phi(), m_vtx0Daug1MassHypo);
360 muM = p4_mu1.M();
361 mupt = p4_mu1.Pt();
362 mueta = p4_mu1.Eta();
363 }
364 MuMass_decor(*mainVertex) = muM;
365 MuPt_decor(*mainVertex) = mupt;
366 MuEta_decor(*mainVertex) = mueta;
367 ChargeMu_decor(*mainVertex) = selectedMuons.at(0)->charge();
368 MuChi2_decor(*mainVertex) = cascadeVertices[1]->trackParticle(0)->chiSquared();
369 MunDoF_decor(*mainVertex) = cascadeVertices[1]->trackParticle(0)->numberDoF();
370
371 // track contribution to chi2 of cascasde fit
372 std::vector< Trk::VxTrackAtVertex > trkAtB = cascadeVertices[1]->vxTrackAtVertex();
373 MuChi2B_decor(*mainVertex) = trkAtB.at(0).trackQuality().chiSquared();
374 MunDoFB_decor(*mainVertex) = trkAtB.at(0).trackQuality().numberDoF();
375
376 //--------------------- Dx / Lambda_c
377 //tagDp = true by default, for Ds+/-, D+ and Lambda_c+
378 float massKX1 = 0.;
379 float massKX1X2 = 0.;
380 float RapidityKX1X2 = 0.;
381 if (dxVertex) {
382 TLorentzVector p4_h1, p4_h2, p4_h3;
383 if(tagDp){
384 p4_h1.SetPtEtaPhiM(dxVertex->trackParticle(0)->pt(),
385 dxVertex->trackParticle(0)->eta(),
386 dxVertex->trackParticle(0)->phi(), m_vtx1Daug1MassHypo);
387 ChargeX1_decor(*mainVertex) = dxVertex->trackParticle(0)->charge();
388 p4_h2.SetPtEtaPhiM(dxVertex->trackParticle(1)->pt(),
389 dxVertex->trackParticle(1)->eta(),
390 dxVertex->trackParticle(1)->phi(), m_vtx1Daug2MassHypo);
391 ChargeK_decor(*mainVertex) = dxVertex->trackParticle(1)->charge();
392 }else{ // Change the order for D- and Lambda_c-
393 p4_h1.SetPtEtaPhiM(dxVertex->trackParticle(0)->pt(),
394 dxVertex->trackParticle(0)->eta(),
395 dxVertex->trackParticle(0)->phi(), m_vtx1Daug2MassHypo);
396 ChargeK_decor(*mainVertex) = dxVertex->trackParticle(0)->charge();
397
398 p4_h2.SetPtEtaPhiM(dxVertex->trackParticle(1)->pt(),
399 dxVertex->trackParticle(1)->eta(),
400 dxVertex->trackParticle(1)->phi(), m_vtx1Daug1MassHypo);
401 ChargeX1_decor(*mainVertex) = dxVertex->trackParticle(1)->charge(); // K (Ds+) or π (D+)
402 }
403 p4_h3.SetPtEtaPhiM(dxVertex->trackParticle(2)->pt(),
404 dxVertex->trackParticle(2)->eta(),
405 dxVertex->trackParticle(2)->phi(), m_vtx1Daug3MassHypo);
406 ChargeX2_decor(*mainVertex) = dxVertex->trackParticle(2)->charge(); // π (Ds+, D+) or p (Lambda_c+)
407
408 massKX1 = (p4_h1 + p4_h2).M();
409 massKX1X2 = (p4_h1 + p4_h2 + p4_h3).M();
410 RapidityKX1X2 = (p4_h1 + p4_h2 + p4_h3).Rapidity();
411 }
412 massKX1_svdecor(*mainVertex) = massKX1;
413 massKX1X2_svdecor(*mainVertex) = massKX1X2;
414 RapidityKX1X2_svdecor(*mainVertex) = RapidityKX1X2;
415
416
417 ATH_CHECK(helper.FillCandwithRefittedVertices(m_refitPV, pvContainer,
418 refPvContainer, &(*m_pvRefitter), m_PV_max, m_DoVertexType, x, 1, mass_b, vtx));
419
420
421 // 4) decorate the main vertex with V0 vertex mass, pt, lifetime and lxy values (plus errors)
422 // V0 points to the main vertex, so lifetime and lxy are w.r.t the main vertex
423 Mass_svdecor(*mainVertex) = m_CascadeTools->invariantMass(moms[0]);
424 MassErr_svdecor(*mainVertex) = m_CascadeTools->invariantMassError(moms[0],x->getCovariance()[0]);
425 Pt_svdecor(*mainVertex) = m_CascadeTools->pT(moms[0]);
426 PtErr_svdecor(*mainVertex) = m_CascadeTools->pTError(moms[0],x->getCovariance()[0]);
427 Lxy_svdecor(*mainVertex) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
428 LxyErr_svdecor(*mainVertex) = m_CascadeTools->lxyError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
429 Tau_svdecor(*mainVertex) = m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1]);
430 TauErr_svdecor(*mainVertex) = m_CascadeTools->tauError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
431
432 // Some checks in DEBUG mode
433 ATH_MSG_DEBUG("chi2 " << x->fitChi2()
434 << " chi2_1 " << m_V0Tools->chisq(cascadeVertices[0])
435 << " chi2_2 " << m_V0Tools->chisq(cascadeVertices[1])
436 << " vprob " << m_CascadeTools->vertexProbability(x->nDoF(),x->fitChi2()));
437 ATH_MSG_DEBUG("ndf " << x->nDoF() << " ndf_1 " << m_V0Tools->ndof(cascadeVertices[0]) << " ndf_2 " << m_V0Tools->ndof(cascadeVertices[1]));
438 ATH_MSG_DEBUG("V0Tools mass_d " << m_V0Tools->invariantMass(cascadeVertices[0],massesDx)
439 << " error " << m_V0Tools->invariantMassError(cascadeVertices[0],massesDx)
440 << " mass_J " << m_V0Tools->invariantMass(cascadeVertices[1],massesMu)
441 << " error " << m_V0Tools->invariantMassError(cascadeVertices[1],massesMu));
442 // masses and errors, using track masses assigned in the fit
443 double Mass_B = m_CascadeTools->invariantMass(moms[1]);
444 double Mass_D = m_CascadeTools->invariantMass(moms[0]);
445 double Mass_B_err = m_CascadeTools->invariantMassError(moms[1],x->getCovariance()[1]);
446 double Mass_D_err = m_CascadeTools->invariantMassError(moms[0],x->getCovariance()[0]);
447 ATH_MSG_DEBUG("Mass_B " << Mass_B << " Mass_D " << Mass_D);
448 ATH_MSG_DEBUG("Mass_B_err " << Mass_B_err << " Mass_D_err " << Mass_D_err);
449 double mprob_B = m_CascadeTools->massProbability(mass_b,Mass_B,Mass_B_err);
450 double mprob_D = m_CascadeTools->massProbability(mass_d,Mass_D,Mass_D_err);
451 ATH_MSG_DEBUG("mprob_B " << mprob_B << " mprob_D " << mprob_D);
452 // masses and errors, assigning user defined track masses
453 ATH_MSG_DEBUG("Mass_b " << m_CascadeTools->invariantMass(moms[1],Masses)
454 << " Mass_d " << m_CascadeTools->invariantMass(moms[0],massesDx));
455 ATH_MSG_DEBUG("Mass_b_err " << m_CascadeTools->invariantMassError(moms[1],x->getCovariance()[1],Masses)
456 << " Mass_d_err " << m_CascadeTools->invariantMassError(moms[0],x->getCovariance()[0],massesDx));
457 ATH_MSG_DEBUG("pt_b " << m_CascadeTools->pT(moms[1])
458 << " pt_d " << m_CascadeTools->pT(moms[0])
459 << " pt_dp " << m_V0Tools->pT(cascadeVertices[0]));
460 ATH_MSG_DEBUG("ptErr_b " << m_CascadeTools->pTError(moms[1],x->getCovariance()[1])
461 << " ptErr_d " << m_CascadeTools->pTError(moms[0],x->getCovariance()[0])
462 << " ptErr_dp " << m_V0Tools->pTError(cascadeVertices[0]));
463 ATH_MSG_DEBUG("lxy_B " << m_V0Tools->lxy(cascadeVertices[1],primaryVertex) << " lxy_D " << m_V0Tools->lxy(cascadeVertices[0],cascadeVertices[1]));
464 ATH_MSG_DEBUG("lxy_b " << m_CascadeTools->lxy(moms[1],cascadeVertices[1],primaryVertex) << " lxy_d " << m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]));
465 ATH_MSG_DEBUG("lxyErr_b " << m_CascadeTools->lxyError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
466 << " lxyErr_d " << m_CascadeTools->lxyError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
467 << " lxyErr_dp " << m_V0Tools->lxyError(cascadeVertices[0],cascadeVertices[1]));
468 ATH_MSG_DEBUG("tau_B " << m_CascadeTools->tau(moms[1],cascadeVertices[1],primaryVertex,mass_b)
469 << " tau_dp " << m_V0Tools->tau(cascadeVertices[0],cascadeVertices[1],massesDx));
470 ATH_MSG_DEBUG("tau_b " << m_CascadeTools->tau(moms[1],cascadeVertices[1],primaryVertex)
471 << " tau_d " << m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1])
472 << " tau_D " << m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1],mass_d));
473 ATH_MSG_DEBUG("tauErr_b " << m_CascadeTools->tauError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
474 << " tauErr_d " << m_CascadeTools->tauError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
475 << " tauErr_dp " << m_V0Tools->tauError(cascadeVertices[0],cascadeVertices[1],massesDx));
476 ATH_MSG_DEBUG("TauErr_b " << m_CascadeTools->tauError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex,mass_b)
477 << " TauErr_d " << m_CascadeTools->tauError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1],mass_d)
478 << " TauErr_dp " << m_V0Tools->tauError(cascadeVertices[0],cascadeVertices[1],massesDx,mass_d));
479
480 ATH_MSG_DEBUG("CascadeTools main vert wrt PV " << " CascadeTools SV " << " V0Tools SV");
481 ATH_MSG_DEBUG("a0z " << m_CascadeTools->a0z(moms[1],cascadeVertices[1],primaryVertex)
482 << ", " << m_CascadeTools->a0z(moms[0],cascadeVertices[0],cascadeVertices[1])
483 << ", " << m_V0Tools->a0z(cascadeVertices[0],cascadeVertices[1]));
484 ATH_MSG_DEBUG("a0zErr " << m_CascadeTools->a0zError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
485 << ", " << m_CascadeTools->a0zError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
486 << ", " << m_V0Tools->a0zError(cascadeVertices[0],cascadeVertices[1]));
487 ATH_MSG_DEBUG("a0xy " << m_CascadeTools->a0xy(moms[1],cascadeVertices[1],primaryVertex)
488 << ", " << m_CascadeTools->a0xy(moms[0],cascadeVertices[0],cascadeVertices[1])
489 << ", " << m_V0Tools->a0xy(cascadeVertices[0],cascadeVertices[1]));
490 ATH_MSG_DEBUG("a0xyErr " << m_CascadeTools->a0xyError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
491 << ", " << m_CascadeTools->a0xyError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
492 << ", " << m_V0Tools->a0xyError(cascadeVertices[0],cascadeVertices[1]));
493 ATH_MSG_DEBUG("a0 " << m_CascadeTools->a0(moms[1],cascadeVertices[1],primaryVertex)
494 << ", " << m_CascadeTools->a0(moms[0],cascadeVertices[0],cascadeVertices[1])
495 << ", " << m_V0Tools->a0(cascadeVertices[0],cascadeVertices[1]));
496 ATH_MSG_DEBUG("a0Err " << m_CascadeTools->a0Error(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
497 << ", " << m_CascadeTools->a0Error(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
498 << ", " << m_V0Tools->a0Error(cascadeVertices[0],cascadeVertices[1]));
499 ATH_MSG_DEBUG("x0 " << m_V0Tools->vtx(cascadeVertices[0]).x() << " y0 " << m_V0Tools->vtx(cascadeVertices[0]).y() << " z0 " << m_V0Tools->vtx(cascadeVertices[0]).z());
500 ATH_MSG_DEBUG("x1 " << m_V0Tools->vtx(cascadeVertices[1]).x() << " y1 " << m_V0Tools->vtx(cascadeVertices[1]).y() << " z1 " << m_V0Tools->vtx(cascadeVertices[1]).z());
501 ATH_MSG_DEBUG("X0 " << primaryVertex->x() << " Y0 " << primaryVertex->y() << " Z0 " << primaryVertex->z());
502 ATH_MSG_DEBUG("rxy0 " << m_V0Tools->rxy(cascadeVertices[0]) << " rxyErr0 " << m_V0Tools->rxyError(cascadeVertices[0]));
503 ATH_MSG_DEBUG("rxy1 " << m_V0Tools->rxy(cascadeVertices[1]) << " rxyErr1 " << m_V0Tools->rxyError(cascadeVertices[1]));
504 ATH_MSG_DEBUG("Rxy0 wrt PV " << m_V0Tools->rxy(cascadeVertices[0],primaryVertex) << " RxyErr0 wrt PV " << m_V0Tools->rxyError(cascadeVertices[0],primaryVertex));
505 ATH_MSG_DEBUG("Rxy1 wrt PV " << m_V0Tools->rxy(cascadeVertices[1],primaryVertex) << " RxyErr1 wrt PV " << m_V0Tools->rxyError(cascadeVertices[1],primaryVertex));
506 ATH_MSG_DEBUG("number of covariance matrices " << (x->getCovariance()).size());
507
508 } // loop over cascadeinfoContainer
509
510 // Deleting cascadeinfo since this won't be stored.
511 // Vertices have been kept in m_cascadeOutputs and should be owned by their container
512 for (auto x : cascadeinfoContainer) delete x;
513
514 return StatusCode::SUCCESS;
515 }
516
517
518 MuPlusDsCascade::MuPlusDsCascade(const std::string& t, const std::string& n, const IInterface* p) : base_class(t,n,p),
519 m_vertexDxContainerKey(""),
520 m_cascadeOutputsKeys{ "MuPlusDsCascadeVtx1", "MuPlusDsCascadeVtx2" },
521 m_VxPrimaryCandidateName("PrimaryVertices"),
522 m_DxMassLower(0.0),
523 m_DxMassUpper(10000.0),
524 m_MassLower(0.0),
525 m_MassUpper(20000.0),
526 m_vtx0MassHypo(-1),
527 m_vtx1MassHypo(-1),
528 m_vtx0Daug1MassHypo(-1),
529 m_vtx1Daug1MassHypo(-1),
530 m_vtx1Daug2MassHypo(-1),
531 m_vtx1Daug3MassHypo(-1),
532 m_Dx_pid(431),
533 m_constrDx(true),
534 m_chi2cut(-1.0),
535 m_iVertexFitter("Trk::TrkVKalVrtFitter"),
536 m_pvRefitter("Analysis::PrimaryVertexRefitter",this),
537 m_V0Tools("Trk::V0Tools"),
538 m_CascadeTools("DerivationFramework::CascadeTools"),
539 m_muonCollectionKey("StacoMuonCollection"),
540 m_trkSelector("InDet::TrackSelectorTool"),
541 m_thresholdPt(0.0),
542 m_mcpCuts(true),
543 m_combOnly(false),
544 m_useCombMeasurement(false)
545 {
546 declareProperty("DxVertices", m_vertexDxContainerKey);
547 declareProperty("VxPrimaryCandidateName", m_VxPrimaryCandidateName);
548 declareProperty("RefPVContainerName", m_refPVContainerName = "RefittedPrimaryVertices");
549 declareProperty("DxMassLowerCut", m_DxMassLower);
550 declareProperty("DxMassUpperCut", m_DxMassUpper);
551 declareProperty("MassLowerCut", m_MassLower);
552 declareProperty("MassUpperCut", m_MassUpper);
553 declareProperty("HypothesisName", m_hypoName = "B");
554 declareProperty("Vtx0MassHypo", m_vtx0MassHypo);
555 declareProperty("Vtx1MassHypo", m_vtx1MassHypo);
556 declareProperty("Vtx0Daug1MassHypo", m_vtx0Daug1MassHypo);
557 declareProperty("Vtx1Daug1MassHypo", m_vtx1Daug1MassHypo);
558 declareProperty("Vtx1Daug2MassHypo", m_vtx1Daug2MassHypo);
559 declareProperty("Vtx1Daug3MassHypo", m_vtx1Daug3MassHypo);
560 declareProperty("DxHypothesis", m_Dx_pid);
561 declareProperty("ApplyDxMassConstraint", m_constrDx);
562 declareProperty("Chi2Cut", m_chi2cut);
563 declareProperty("RefitPV", m_refitPV = true);
564 declareProperty("MaxnPV", m_PV_max = 999);
565 declareProperty("MinNTracksInPV", m_PV_minNTracks = 0);
566 declareProperty("DoVertexType", m_DoVertexType = 7);
567 declareProperty("TrkVertexFitterTool", m_iVertexFitter);
568 declareProperty("PVRefitter", m_pvRefitter);
569 declareProperty("V0Tools", m_V0Tools);
570 declareProperty("CascadeTools", m_CascadeTools);
571 declareProperty("CascadeVertexCollections", m_cascadeOutputsKeys);
572 declareProperty("muonCollectionKey", m_muonCollectionKey);
573 declareProperty("TrackSelectorTool", m_trkSelector);
574 declareProperty("muonThresholdPt", m_thresholdPt);
575 declareProperty("useMCPCuts", m_mcpCuts);
576 declareProperty("combOnly", m_combOnly);
577 declareProperty("useCombinedMeasurement", m_useCombMeasurement);
578 }
579
580 MuPlusDsCascade::~MuPlusDsCascade(){ }
581
582 StatusCode MuPlusDsCascade::performSearch(std::vector<Trk::VxCascadeInfo*> *cascadeinfoContainer, const EventContext&) const
583 {
584 ATH_MSG_DEBUG( "MuPlusDsCascade::performSearch" );
585 assert(cascadeinfoContainer!=nullptr);
586
587 // Get TrackParticle container (for setting links to the original tracks)
588 const xAOD::TrackParticleContainer *trackContainer{};
589 ATH_CHECK(evtStore()->retrieve(trackContainer , "InDetTrackParticles" )); // FIXME Use Handles
590
591 // Get V0 container
592 const xAOD::VertexContainer *dxContainer{};
593 ATH_CHECK(evtStore()->retrieve(dxContainer , m_vertexDxContainerKey )); // FIXME Use Handles
594
595 double mass_d = m_vtx1MassHypo;
596 std::vector<const xAOD::TrackParticle*> tracksMu;
597 std::vector<const xAOD::TrackParticle*> tracksDx;
598 std::vector<double> massesMu;
599 massesMu.push_back(m_vtx0Daug1MassHypo); // µ
600 std::vector<double> massesDx;
601 massesDx.push_back(m_vtx1Daug1MassHypo); // K (Ds+/-) or π (D+, Lambda_c+)
602 massesDx.push_back(m_vtx1Daug2MassHypo); // K
603 massesDx.push_back(m_vtx1Daug3MassHypo); // π (Ds+/-, D+) ot p (Lambda_c+)
604 std::vector<double> massesDm; // Alter the order of masses for D- and Lambda_c-
605 massesDm.push_back(m_vtx1Daug2MassHypo); // K
606 massesDm.push_back(m_vtx1Daug1MassHypo); // π
607 massesDm.push_back(m_vtx1Daug3MassHypo); // π (D-) or p (Lambda_c-)
608 std::vector<double> Masses;
609 Masses.push_back(m_vtx0Daug1MassHypo); // µ
610 Masses.push_back(m_vtx1MassHypo); // Dx / Lambda_c
611
612 //-------------------------------------------------------------------------------
613 // Retrieving and selecting muons (taken from JpsiFinder.cxx)
614 // Get the muons from StoreGate
615 const xAOD::MuonContainer* importedMuonCollection{};
616 StatusCode sc = evtStore()->retrieve(importedMuonCollection,m_muonCollectionKey); // FIXME Use Handles
617
618 if(sc.isFailure()){
619 ATH_MSG_WARNING("No muon collection with key " << m_muonCollectionKey << " found in StoreGate");
620 return StatusCode::SUCCESS;;
621 }else{
622 ATH_MSG_DEBUG("Found muon collections with key "<<m_muonCollectionKey);
623 }
624 ATH_MSG_DEBUG("Muon container size "<<importedMuonCollection->size());
625
626 // Typedef for vectors of muons
627 typedef std::vector<const xAOD::Muon*> MuonBag;
628
629 // Select the muons
630 const xAOD::Vertex* vx = 0;
631 MuonBag theMuonsAfterSelection;
632 for (auto mu : *importedMuonCollection) {
633 if ( !mu ) continue;
634 if (!mu->inDetTrackParticleLink().isValid()) continue; // No muons without ID tracks
635 const xAOD::TrackParticle* muonTrk = *(mu->inDetTrackParticleLink());
636 if ( !muonTrk) continue;
637 if ( !m_trkSelector->decision(*muonTrk, vx) ) continue; // all ID tracks must pass basic tracking cuts
638 if ( std::fabs(muonTrk->pt())<m_thresholdPt ) continue; // higher pt cut if needed
639 if ( m_mcpCuts && !mu->passesIDCuts()) continue; // cuts of the MCP group recommendation
640 if ( m_combOnly && mu->muonType() != xAOD::Muon::Combined ) continue; // require combined muons
641 if ( mu->muonType() == xAOD::Muon::SiliconAssociatedForwardMuon && !m_useCombMeasurement) continue;
642
643 theMuonsAfterSelection.push_back(mu);
644 }
645 if (theMuonsAfterSelection.size() == 0) return StatusCode::SUCCESS;;
646 ATH_MSG_DEBUG("Number of muons after selection: " << theMuonsAfterSelection.size());
647 //-------------------------------------------------------------------------------
648
649 //-------------------------------------------------------------------------------
650 // Select the D_s+/D+/Lambda_c+ candidates before calling cascade fit
651 std::vector<const xAOD::Vertex*> selectedDxCandidates;
652
653
654 for(auto vxcItr : *dxContainer){
655
656 // Check the passed flag first
657 const xAOD::Vertex* vtx = vxcItr;
658 if(std::abs(m_Dx_pid)==431) { // D_s+/-
659 SG::AuxElement::Accessor<Char_t> flagAcc1("passed_Ds");
660 if(flagAcc1.isAvailable(*vtx)){
661 if(!flagAcc1(*vtx)) continue;
662 }
663 }
664
665 if(std::abs(m_Dx_pid)==411) { // D+/-
666 SG::AuxElement::Accessor<Char_t> flagAcc1("passed_Dp");
667 SG::AuxElement::Accessor<Char_t> flagAcc2("passed_Dm");
668 bool isDp(true);
669 bool isDm(true);
670 if(flagAcc1.isAvailable(*vtx)){
671 if(!flagAcc1(*vtx)) isDp = false;
672 }
673 if(flagAcc2.isAvailable(*vtx)){
674 if(!flagAcc2(*vtx)) isDm = false;
675 }
676 if(!(isDp||isDm)) continue;
677 }
678
679 // Track selection - Loose
680 if ( !m_trackSelectionTools->accept(vxcItr->trackParticle(0)) ){
681 ATH_MSG_DEBUG(" Original Dx/Lambda_c candidate rejected by the track's cut level - loose ");
682 continue;
683 }
684 if ( !m_trackSelectionTools->accept(vxcItr->trackParticle(1)) ){
685 ATH_MSG_DEBUG(" Original Dx/Lambda_c candidate rejected by the track's cut level - loose ");
686 continue;
687 }
688 if ( !m_trackSelectionTools->accept(vxcItr->trackParticle(2)) ){
689 ATH_MSG_DEBUG(" Original Dx/Lambda_c candidate rejected by the track's cut level - loose ");
690 continue;
691 }
692
693 // Ensure the total charge is correct
694 if( std::abs( vxcItr->trackParticle(0)->charge()+vxcItr->trackParticle(1)->charge()+vxcItr->trackParticle(2)->charge() ) != 1 ){
695 ATH_MSG_DEBUG(" Original Dx/Lambda_c candidate rejected by the charge requirement: "
696 << vxcItr->trackParticle(0)->charge() << ", " << vxcItr->trackParticle(1)->charge() << ", " << vxcItr->trackParticle(2)->charge() );
697 continue;
698 }
699
700 // Check D_(s)/Lambda_c +/- candidate invariant mass and skip if need be
701 double mass_D;
702 if( (std::abs(m_Dx_pid) == 411 || std::abs(m_Dx_pid) == 4122) && vxcItr->trackParticle(2)->charge()<0) // D- & Lambda_c-
703 mass_D = m_V0Tools->invariantMass(vxcItr,massesDm);
704 else // D+, D_s+/-, Lambda_c+
705 mass_D = m_V0Tools->invariantMass(vxcItr,massesDx);
706 ATH_MSG_DEBUG("D_(s)/Lambda_c mass " << mass_D);
707 if(mass_D < m_DxMassLower || mass_D > m_DxMassUpper) {
708 ATH_MSG_DEBUG(" Original D_(s)/Lambda_c candidate rejected by the mass cut: mass = "
709 << mass_D << " != (" << m_DxMassLower << ", " << m_DxMassUpper << ")" );
710 continue;
711 }
712
713 // Add loose cut on phi(K+K-) mass for D_s->phi π
714 if(std::abs(m_Dx_pid)==431){
715 TLorentzVector p4Kp_in, p4Km_in;
716 p4Kp_in.SetPtEtaPhiM( vxcItr->trackParticle(0)->pt(),
717 vxcItr->trackParticle(0)->eta(),
718 vxcItr->trackParticle(0)->phi(), m_vtx1Daug1MassHypo);
719 p4Km_in.SetPtEtaPhiM( vxcItr->trackParticle(1)->pt(),
720 vxcItr->trackParticle(1)->eta(),
721 vxcItr->trackParticle(1)->phi(), m_vtx1Daug2MassHypo);
722 double mass_phi = (p4Kp_in + p4Km_in).M();
723 ATH_MSG_DEBUG("phi mass " << mass_phi);
724 if(mass_phi > 1600) { //1st cut on phi(K+K-) mass, loose one
725 ATH_MSG_DEBUG(" Original phi candidate rejected by the mass cut: mass = " << mass_phi );
726 continue;
727
728 }
729 }
730
731 // Cut on pT(p) for Lambda_c
732 if(std::abs(m_Dx_pid)==4122){
733 if (vxcItr->trackParticle(2)->pt() < 2400) continue;
734 }
735
736 // L.G., 26.03.21: Check there are no previous D+/- combinations with inverse pion positions
737 if(std::abs(m_Dx_pid)==411 && selectedDxCandidates.size()>0) {
738 bool Dpmcopy(false);
739 for(auto dxItr : selectedDxCandidates){
740 if(vxcItr->trackParticle(2)/*π*/->charge()>0) { // D+
741 if ( vxcItr->trackParticle(2) == dxItr->trackParticle(0) && vxcItr->trackParticle(0) == dxItr->trackParticle(2)
742 && vxcItr->trackParticle(1) == dxItr->trackParticle(1) ) Dpmcopy = true;
743 } else { // D-
744 if ( vxcItr->trackParticle(2) == dxItr->trackParticle(1) && vxcItr->trackParticle(1) == dxItr->trackParticle(2)
745 && vxcItr->trackParticle(0) == dxItr->trackParticle(0) ) Dpmcopy = true;
746 }
747 }
748 if (Dpmcopy) continue;
749 }
750
751 selectedDxCandidates.push_back(vxcItr);
752 } // end for(auto vxcItr : dxContainer)
753 if(selectedDxCandidates.size()<1) return StatusCode::SUCCESS;
754 //-------------------------------------------------------------------------------
755
756 //-------------------------------------------------------------------------------
757 // Select mu D_(s)+/Lambda_c candidates
758 // Iterate over muons
759 for(auto muItr : theMuonsAfterSelection){
760 tracksMu.clear();
761 //Convert to trackParticle base
762 auto TrkMuon = muItr->trackParticle( xAOD::Muon::InnerDetectorTrackParticle );
763 tracksMu.push_back(TrkMuon);
764 if (tracksMu.size() != 1 || massesMu.size() != 1 ) {
765 ATH_MSG_WARNING("Problems with muon input");
766 }
767
768 // Iterate over D_(s)/Lambda_c +/- vertices
769 for(auto dxItr : selectedDxCandidates){
770 // Check identical tracks in input
771 if(std::find(tracksMu.cbegin(), tracksMu.cend(), dxItr->trackParticle(0)) != tracksMu.cend()) continue;
772 if(std::find(tracksMu.cbegin(), tracksMu.cend(), dxItr->trackParticle(1)) != tracksMu.cend()) continue;
773 if(std::find(tracksMu.cbegin(), tracksMu.cend(), dxItr->trackParticle(2)) != tracksMu.cend()) continue;
774
775 size_t dxTrkNum = dxItr->nTrackParticles();
776 tracksDx.clear();
777 for( unsigned int it=0; it<dxTrkNum; it++) tracksDx.push_back(dxItr->trackParticle(it));
778 if (tracksDx.size() != 3 || massesDx.size() != 3 ) {
779 ATH_MSG_WARNING("Problems with D_(s)/Lambda_c +/- input");
780 }
781
782 if(std::find(tracksDx.cbegin(), tracksDx.cend(), TrkMuon) != tracksDx.cend()) continue; //looking for the muons in the D+ tracks
783
784 ATH_MSG_DEBUG("Using tracks" << tracksMu[0] << ", " << tracksDx[0] << ", " << tracksDx[1] << ", " << tracksDx[2]);
785 // Apply the user's settings to the fitter
786 // Reset
787 std::unique_ptr<Trk::IVKalState> state (m_iVertexFitter->makeState());
788 // Robustness
789 int robustness = 0;
790 m_iVertexFitter->setRobustness(robustness, *state);
791 // Build up the topology
792 // Vertex list
793 std::vector<Trk::VertexID> vrtList;
794 // D_(s)/Lambda_c +/- vertex
795 Trk::VertexID vID;
796 if (m_constrDx) {
797 if((std::abs(m_Dx_pid)==411 || std::abs(m_Dx_pid)==4122 ) && dxItr->trackParticle(2)->charge()<0) // D-, Lambda_c-
798 vID = m_iVertexFitter->startVertex(tracksDx,massesDm, *state,mass_d);
799 else // D+, D_s+/-, Lambda_c+
800 vID = m_iVertexFitter->startVertex(tracksDx,massesDx, *state,mass_d);
801 } else {
802 if((std::abs(m_Dx_pid)==411 || std::abs(m_Dx_pid)==4122 ) && dxItr->trackParticle(2)->charge()<0) // D-, Lambda_c-
803 vID = m_iVertexFitter->startVertex(tracksDx,massesDm, *state);
804 else // D+, D_s+/-, Lambda_c+
805 vID = m_iVertexFitter->startVertex(tracksDx,massesDx, *state);
806 }
807 vrtList.push_back(vID);
808
809 // B vertex including muon
810 //For mass constraint, not used yet
811 Trk::VertexID vID2 = m_iVertexFitter->nextVertex(tracksMu,massesMu,vrtList, *state);
812 ATH_MSG_DEBUG(vID2);
813
814 // Do the work
815 std::unique_ptr<Trk::VxCascadeInfo> result(m_iVertexFitter->fitCascade(*state));
816
817 if (result) {
818 // reset links to original tracks
819 BPhysPVCascadeTools::PrepareVertexLinks(result.get(), trackContainer);
820 ATH_MSG_DEBUG("storing tracks " << ((result->vertices())[0] /*D_(s)+*/)->trackParticle(0) << ", "
821 << ((result->vertices())[0])->trackParticle(1) << ", "
822 << ((result->vertices())[0])->trackParticle(2) << ", "
823 << ((result->vertices())[1]/*mu*/)->trackParticle(0));
824
825 // necessary to prevent memory leak
826 result->setSVOwnership(true);
827
828 // Chi2/DOF cut
829 double bChi2DOF = result->fitChi2()/result->nDoF();
830 ATH_MSG_DEBUG("Candidate chi2/DOF is " << bChi2DOF);
831 bool chi2CutPassed = (m_chi2cut <= 0.0 || bChi2DOF < m_chi2cut);
832
833 const std::vector<xAOD::Vertex*> &cascadeVertices = result->vertices();
834 const std::vector< std::vector<TLorentzVector> > &moms = result->getParticleMoms();
835
836 //----------------------------------------------------
837 // retrieve primary vertices
838 //----------------------------------------------------
839 const xAOD::Vertex * primaryVertex{};
840 const xAOD::VertexContainer *pvContainer{};
841 ATH_CHECK(evtStore()->retrieve(pvContainer, m_VxPrimaryCandidateName)); // FIXME Use Handles
842 ATH_MSG_DEBUG("Found " << m_VxPrimaryCandidateName << " in StoreGate!");
843
844 if (pvContainer->size()==0){
845 ATH_MSG_WARNING("You have no primary vertices: " << pvContainer->size());
846 return StatusCode::RECOVERABLE;
847 } else {
848 primaryVertex = (*pvContainer)[0];
849 }
850
851 // Add stronger cut on phi(K+K-) mass for D_s->phi π after cascade fit
852 if(std::abs(m_Dx_pid)==431){
853 TLorentzVector p4Kp_in, p4Km_in;
854 p4Kp_in.SetPtEtaPhiM( cascadeVertices[0]->trackParticle(0)->pt(),
855 cascadeVertices[0]->trackParticle(0)->eta(),
856 cascadeVertices[0]->trackParticle(0)->phi(), m_vtx1Daug1MassHypo);
857 p4Km_in.SetPtEtaPhiM( cascadeVertices[0]->trackParticle(1)->pt(),
858 cascadeVertices[0]->trackParticle(1)->eta(),
859 cascadeVertices[0]->trackParticle(1)->phi(), m_vtx1Daug2MassHypo);
860 double mass_phi = (p4Kp_in + p4Km_in).M();
861 ATH_MSG_DEBUG("phi mass " << mass_phi);
862 if(mass_phi > 1100) {
863 ATH_MSG_DEBUG(" Original phi candidate rejected by the mass cut: mass = " << mass_phi );
864 continue;
865 }
866 }
867
868 double mass = m_CascadeTools->invariantMass(moms[1]);
869 if(chi2CutPassed) {
870 if (mass >= m_MassLower && mass <= m_MassUpper) {
871 if (m_CascadeTools->lxy(moms[1],cascadeVertices[1],primaryVertex) > 0.09){ //B_Lxy
872 if (m_CascadeTools->pT(moms[1]) > 9500){ //B_pT
873 if (m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]) > 0){ //Dx_lxy
874 if (std::abs(m_Dx_pid)!=411) cascadeinfoContainer->push_back(result.release());
875 else if (m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]) > 0.09){ //D+_Lxy
876 cascadeinfoContainer->push_back(result.release());
877 }//D+_Lxy
878 }//Dx_lxy
879 } //B_pT
880 }//B_Lxy
881 } else {
882 ATH_MSG_DEBUG("Candidate rejected by the mass cut: mass = "
883 << mass << " != (" << m_MassLower << ", " << m_MassUpper << ")" );
884 }
885 } // chi2CutPassed
886
887 } // result != nullptr
888
889 } //Iterate over D_(s)/Lambda_c +/- vertices
890
891 } //Iterate over muons
892
893 ATH_MSG_DEBUG("cascadeinfoContainer size " << cascadeinfoContainer->size());
894
895 return StatusCode::SUCCESS;
896 }
897
898}
899
900
eta
Scalar eta() const
pseudorapidity method
Definition AmgMatrixBasePlugin.h:83
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
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_INFO
#define ATH_MSG_INFO(x)
Definition AthMsgStreamMacros.h:31
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition AthMsgStreamMacros.h:32
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
BPHYS_CHECK
#define BPHYS_CHECK(EXP)
Useful CHECK macro.
Definition BPhysHelper.h:738
BPhysHypoHelper.h
: B-physics xAOD helpers.
ConstDataVector.h
DataVector adapter that acts like it holds const pointers.
ANA_CHECK
#define ANA_CHECK(EXP)
check whether the given expression was successful
Definition Control/AthToolSupport/AsgMessaging/AsgMessaging/MessageCheck.h:324
HepMCHelpers.h
ATLAS-specific HepMC functions.
sc
static Double_t sc
Definition LArPhysWaveHECTool.cxx:37
muonContainer
xAOD::MuonContainer * muonContainer
Definition TrigGlobEffCorrValidation.cxx:188
x
#define x
DataVector::size
size_type size() const noexcept
Returns the number of elements in the collection.
DerivationFramework::BPhysPVCascadeTools::LinkVertices
static bool LinkVertices(SG::AuxElement::Decorator< VertexLinkVector > &decor, const std::vector< const xAOD::Vertex * > &vertices, const xAOD::VertexContainer *vertexContainer, const xAOD::Vertex *vert)
Definition BPhysPVCascadeTools.cxx:460
DerivationFramework::BPhysPVCascadeTools::PrepareVertexLinks
static void PrepareVertexLinks(Trk::VxCascadeInfo *result, const xAOD::TrackParticleContainer *importedTrackCollection)
Definition BPhysPVCascadeTools.cxx:204
DerivationFramework::BPhysPVCascadeTools::getParticleMass
static double getParticleMass(const HepPDT::ParticleDataTable *pdt, int pdg)
Definition BPhysPVCascadeTools.cxx:491
DerivationFramework::BPhysPVCascadeTools::FindVertex
static const xAOD::Vertex * FindVertex(const xAOD::VertexContainer *c, const xAOD::Vertex *v)
Definition BPhysPVCascadeTools.h:145
DerivationFramework::BPhysPVCascadeTools::SetVectorInfo
static void SetVectorInfo(xAOD::BPhysHelper &, const Trk::VxCascadeInfo *)
Definition BPhysPVCascadeTools.cxx:424
DerivationFramework::MuPlusDsCascade::m_trkSelector
ToolHandle< Trk::ITrackSelectorTool > m_trkSelector
Definition MuPlusDsCascade.h:94
DerivationFramework::MuPlusDsCascade::m_eventInfo_key
SG::ReadHandleKey< xAOD::EventInfo > m_eventInfo_key
Definition MuPlusDsCascade.h:76
DerivationFramework::MuPlusDsCascade::m_VxPrimaryCandidateName
std::string m_VxPrimaryCandidateName
Name of primary vertex container.
Definition MuPlusDsCascade.h:59
DerivationFramework::MuPlusDsCascade::m_cascadeOutputsKeys
std::vector< std::string > m_cascadeOutputsKeys
Definition MuPlusDsCascade.h:57
DerivationFramework::MuPlusDsCascade::m_CascadeTools
ToolHandle< DerivationFramework::CascadeTools > m_CascadeTools
Definition MuPlusDsCascade.h:80
DerivationFramework::MuPlusDsCascade::addBranches
virtual StatusCode addBranches(const EventContext &ctx) const override
Definition MuPlusDsCascade.cxx:95
DerivationFramework::MuPlusDsCascade::m_V0Tools
ToolHandle< Trk::V0Tools > m_V0Tools
Definition MuPlusDsCascade.h:79
DerivationFramework::MuPlusDsCascade::initialize
virtual StatusCode initialize() override
Definition MuPlusDsCascade.cxx:38
DerivationFramework::MuPlusDsCascade::m_pvRefitter
ToolHandle< Analysis::PrimaryVertexRefitter > m_pvRefitter
Definition MuPlusDsCascade.h:78
DerivationFramework::MuPlusDsCascade::m_trackSelectionTools
std::unique_ptr< InDet::InDetTrackSelectionTool > m_trackSelectionTools
Definition MuPlusDsCascade.h:81
DerivationFramework::MuPlusDsCascade::MuPlusDsCascade
MuPlusDsCascade(const std::string &t, const std::string &n, const IInterface *p)
Definition MuPlusDsCascade.cxx:518
DerivationFramework::MuPlusDsCascade::performSearch
StatusCode performSearch(std::vector< Trk::VxCascadeInfo * > *cascadeinfoContainer, const EventContext &ctx) const
Definition MuPlusDsCascade.cxx:582
DerivationFramework::MuPlusDsCascade::m_partPropSvc
ServiceHandle< IPartPropSvc > m_partPropSvc
Definition MuPlusDsCascade.h:82
DerivationFramework::MuPlusDsCascade::m_iVertexFitter
ToolHandle< Trk::TrkVKalVrtFitter > m_iVertexFitter
Definition MuPlusDsCascade.h:77
DerivationFramework::MuPlusDsCascade::m_hypoName
std::string m_hypoName
name of the mass hypothesis.
Definition MuPlusDsCascade.h:86
SG::AuxElement::Decorator
SG::Decorator< T, ALLOC > Decorator
Definition AuxElement.h:575
SG::AuxElement::Accessor
SG::Accessor< T, ALLOC > Accessor
Definition AuxElement.h:572
SG::ConstAccessor::isAvailable
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
xAOD::BPhysHypoHelper::setMass
bool setMass(const float val)
Set given invariant mass and its error.
Definition BPhysHypoHelper.cxx:52
xAOD::BPhysHypoHelper::setMassErr
bool setMassErr(const float val)
invariant mass error
Definition BPhysHypoHelper.cxx:57
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::TrackParticle_v1::charge
float charge() const
Returns the charge.
Definition TrackParticle_v1.cxx:143
xAOD::Vertex_v1::z
float z() const
Returns the z position.
xAOD::Vertex_v1::trackParticle
const TrackParticle * trackParticle(size_t i) const
Get the pointer to a given track that was used in vertex reco.
Definition Vertex_v1.cxx:249
xAOD::Vertex_v1::y
float y() const
Returns the y position.
xAOD::Vertex_v1::x
float x() const
Returns the x position.
contains
bool contains(const std::string &s, const std::string &regx)
does a string contain the substring
Definition hcg.cxx:114
DerivationFramework
THE reconstruction tool.
Definition ParticleSortingAlg.h:24
DerivationFramework::MuonsLink
ElementLink< xAOD::MuonContainer > MuonsLink
Definition MuPlusDsCascade.cxx:35
DerivationFramework::TrackBag
std::vector< const xAOD::TrackParticle * > TrackBag
Definition JpsiPlusDpstCascade.cxx:25
DerivationFramework::VertexLink
ElementLink< xAOD::VertexContainer > VertexLink
Definition JpsiPlusDpstCascade.cxx:23
DerivationFramework::MuonsLinkVector
std::vector< MuonsLink > MuonsLinkVector
Definition MuPlusDsCascade.cxx:36
DerivationFramework::VertexLinkVector
std::vector< VertexLink > VertexLinkVector
Definition JpsiPlusDpstCascade.cxx:24
MC::KPLUS
static const int KPLUS
Definition HepMCHelpers.h:117
MC::MUON
static const int MUON
Definition HepMCHelpers.h:78
MC::DPLUS
static const int DPLUS
Definition HepMCHelpers.h:118
MC::BCPLUS
static const int BCPLUS
Definition HepMCHelpers.h:124
MC::DSPLUS
static const int DSPLUS
Definition HepMCHelpers.h:121
MC::PIPLUS
static const int PIPLUS
Definition HepMCHelpers.h:110
MC::LAMBDACPLUS
static const int LAMBDACPLUS
Definition HepMCHelpers.h:128
MC::PROTON
static const int PROTON
Definition HepMCHelpers.h:125
xAOD::VertexAuxContainer
VertexAuxContainer_v1 VertexAuxContainer
Definition of the current jet auxiliary container.
Definition VertexAuxContainer.h:19
xAOD::TrackParticle
TrackParticle_v1 TrackParticle
Reference the current persistent version:
Definition Event/xAOD/xAODTracking/xAODTracking/TrackParticle.h:13
xAOD::VertexContainer
VertexContainer_v1 VertexContainer
Definition of the current "Vertex container version".
Definition VertexContainer.h:14
xAOD::Vertex
Vertex_v1 Vertex
Define the latest version of the vertex class.
Definition Event/xAOD/xAODTracking/xAODTracking/Vertex.h:16
xAOD::TrackParticleContainer
TrackParticleContainer_v1 TrackParticleContainer
Definition of the current "TrackParticle container version".
Definition Event/xAOD/xAODTracking/xAODTracking/TrackParticleContainer.h:14
xAOD::Muon
Muon_v1 Muon
Reference the current persistent version:
Definition Event/xAOD/xAODMuon/xAODMuon/Muon.h:13
xAOD::MuonContainer
MuonContainer_v1 MuonContainer
Definition of the current "Muon container version".
Definition Event/xAOD/xAODMuon/xAODMuon/MuonContainer.h:14
result
Generated on for ATLAS Offline Software by doxygen 1.14.0