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