ATLAS Offline Software
Loading...
Searching...
No Matches
MuPlusDpstCascade.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3*/
5// MuPlusDpstCascade.cxx, (c) ATLAS Detector software
7#include "MuPlusDpstCascade.h"
8#include "TrkVertexFitterInterfaces/IVertexFitter.h"
9#include "TrkVKalVrtFitter/TrkVKalVrtFitter.h"
10#include "TrkVertexAnalysisUtils/V0Tools.h"
11#include "DerivationFrameworkBPhys/CascadeTools.h"
12#include "BPhysPVCascadeTools.h"
13#include "xAODTracking/VertexAuxContainer.h"
14#include "xAODTracking/Vertex.h"
15#include "xAODBPhys/BPhysHypoHelper.h"
16#include <algorithm>
17#include "xAODTracking/VertexContainer.h"
18#include "GeneratorModules/GenData.h"
19#include "TrkVKalVrtFitter/VxCascadeInfo.h"
20#include "TruthUtils/HepMCHelpers.h"
21
22namespace DerivationFramework {
23 typedef ElementLink<xAOD::VertexContainer> VertexLink;
24 typedef std::vector<VertexLink> VertexLinkVector;
25 typedef std::vector<const xAOD::TrackParticle*> TrackBag;
26
27
28 StatusCode MuPlusDpstCascade::initialize() {
29 // retrieving vertex Fitter
30 ATH_CHECK( m_iVertexFitter.retrieve());
31
32 // retrieving the V0 tools
33 ATH_CHECK( m_V0Tools.retrieve());
34
35 // retrieving the Cascade tools
36 ATH_CHECK( m_CascadeTools.retrieve());
37
38 ATH_CHECK( m_eventInfo_key.initialize() );
39
40 //======================== inDetTrack selection tool ==================
41 m_trackSelectionTools = std::make_unique<InDet::InDetTrackSelectionTool>("TrackSelector");
42 ANA_CHECK(m_trackSelectionTools->setProperty("CutLevel", "LoosePrimary"));
43 ANA_CHECK(m_trackSelectionTools->initialize() );
44 //=====================================================================
45
46 auto gendata = std::make_shared<GenData>();
47 if(m_vtx0MassHypo < 0.) m_vtx0MassHypo = gendata->particleMass( MC::BCPLUS).value();
48 if(m_vtx1MassHypo < 0.) m_vtx1MassHypo = gendata->particleMass( MC::D0).value();
49 if(m_vtx0Daug1MassHypo < 0.) m_vtx0Daug1MassHypo = gendata->particleMass( MC::MUON).value();
50 if(m_vtx0Daug2MassHypo < 0.) m_vtx0Daug2MassHypo = gendata->particleMass( MC::PIPLUS).value();
51 if(m_vtx1Daug1MassHypo < 0.) m_vtx1Daug1MassHypo = gendata->particleMass( MC::PIPLUS).value();
52 if(m_vtx1Daug2MassHypo < 0.) m_vtx1Daug2MassHypo = gendata->particleMass( MC::KPLUS).value();
53
54 return StatusCode::SUCCESS;
55 }
56
57
58 StatusCode MuPlusDpstCascade::addBranches(const EventContext& ctx) const
59 {
60 std::vector<Trk::VxCascadeInfo*> cascadeinfoContainer;
61 constexpr int topoN = 2;
62 std::array<xAOD::VertexContainer*, topoN> Vtxwritehandles;
63 std::array<xAOD::VertexAuxContainer*, topoN> Vtxwritehandlesaux;
64 if(m_cascadeOutputsKeys.size() !=topoN) { ATH_MSG_FATAL("Incorrect number of VtxContainers"); return StatusCode::FAILURE; }
65
66 for(int i =0; i<topoN;i++){
67 Vtxwritehandles[i] = new xAOD::VertexContainer();
68 Vtxwritehandlesaux[i] = new xAOD::VertexAuxContainer();
69 Vtxwritehandles[i]->setStore(Vtxwritehandlesaux[i]);
70 ATH_CHECK(evtStore()->record(Vtxwritehandles[i] , m_cascadeOutputsKeys[i] )); // FIXME Use Handles
71 ATH_CHECK(evtStore()->record(Vtxwritehandlesaux[i], m_cascadeOutputsKeys[i] + "Aux.")); // FIXME Use Handles
72 }
73
74 //----------------------------------------------------
75 // retrieve primary vertices
76 //----------------------------------------------------
77 const xAOD::Vertex * primaryVertex{};
78 const xAOD::VertexContainer *pvContainer{};
79 ATH_CHECK(evtStore()->retrieve(pvContainer, m_VxPrimaryCandidateName)); // FIXME Use Handles
80 ATH_MSG_DEBUG("Found " << m_VxPrimaryCandidateName << " in StoreGate!");
81
82 if (pvContainer->size()==0){
83 ATH_MSG_WARNING("You have no primary vertices: " << pvContainer->size());
84 return StatusCode::RECOVERABLE;
85 } else {
86 primaryVertex = (*pvContainer)[0];
87 }
88
89 //----------------------------------------------------
90 // Try to retrieve refitted primary vertices
91 //----------------------------------------------------
92 xAOD::VertexContainer* refPvContainer{};
93 xAOD::VertexAuxContainer* refPvAuxContainer{};
94 if (m_refitPV) {
95 if (evtStore()->contains<xAOD::VertexContainer>(m_refPVContainerName)) {
96 // refitted PV container exists. Get it from the store gate
97 ATH_CHECK(evtStore()->retrieve(refPvContainer , m_refPVContainerName )); // FIXME Use Handles
98 ATH_CHECK(evtStore()->retrieve(refPvAuxContainer, m_refPVContainerName + "Aux.")); // FIXME Use Handles
99 } else {
100 // refitted PV container does not exist. Create a new one.
101 refPvContainer = new xAOD::VertexContainer;
102 refPvAuxContainer = new xAOD::VertexAuxContainer;
103 refPvContainer->setStore(refPvAuxContainer);
104 ATH_CHECK(evtStore()->record(refPvContainer , m_refPVContainerName)); // FIXME Use Handles
105 ATH_CHECK(evtStore()->record(refPvAuxContainer, m_refPVContainerName+"Aux.")); // FIXME Use Handles
106 }
107 }
108
109 ATH_CHECK(performSearch(&cascadeinfoContainer, ctx));
110 SG::ReadHandle<xAOD::EventInfo> evt(m_eventInfo_key, ctx);
111 if(!evt.isValid()) {
112 ATH_MSG_ERROR("Cannot Retrieve " << m_eventInfo_key.key() );
113 return StatusCode::FAILURE;
114 }
115
116 BPhysPVCascadeTools helper(&(*m_CascadeTools), evt.cptr());
117 helper.SetMinNTracksInPV(m_PV_minNTracks);
118 // Decorators for the main vertex: chi2, ndf, pt and pt error, plus the D0 vertex variables
119 SG::AuxElement::Decorator<VertexLinkVector> CascadeLinksDecor("CascadeVertexLinks");
120 SG::AuxElement::Decorator<VertexLinkVector> MuPiLinksDecor("MuPiVertexLinks");
121 SG::AuxElement::Decorator<VertexLinkVector> D0LinksDecor("D0VertexLinks");
122 SG::AuxElement::Decorator<float> chi2_decor("ChiSquared");
123 SG::AuxElement::Decorator<float> ndof_decor("NumberDoF");
124 SG::AuxElement::Decorator<float> Pt_decor("Pt");
125 SG::AuxElement::Decorator<float> PtErr_decor("PtErr");
126 SG::AuxElement::Decorator<float> Mass_svdecor("D0_mass");
127 SG::AuxElement::Decorator<float> MassErr_svdecor("D0_massErr");
128 SG::AuxElement::Decorator<float> Pt_svdecor("D0_Pt");
129 SG::AuxElement::Decorator<float> PtErr_svdecor("D0_PtErr");
130 SG::AuxElement::Decorator<float> Lxy_svdecor("D0_Lxy");
131 SG::AuxElement::Decorator<float> LxyErr_svdecor("D0_LxyErr");
132 SG::AuxElement::Decorator<float> Tau_svdecor("D0_Tau");
133 SG::AuxElement::Decorator<float> TauErr_svdecor("D0_TauErr");
134
135 SG::AuxElement::Decorator<float> massMuPi_decor("MuPi_mass"); //mu+pi_soft mass before fit
136 SG::AuxElement::Decorator<float> MassKpi_svdecor("Kpi_mass");
137 SG::AuxElement::Decorator<float> MassMuPiAft_decor("MuPiAft_mass"); //mu+pi_soft mass after fit
138 SG::AuxElement::Decorator<float> MassPiD0_decor("PiD0_mass");
139 SG::AuxElement::Decorator<float> MassMuPiPiK_decor("MuPiPiK_mass");
140
141 SG::AuxElement::Decorator<float> ChargePi_decor("Pi_charge");
142 SG::AuxElement::Decorator<float> ChargeMu_decor("Mu_charge");
143 SG::AuxElement::Decorator<float> ChargeK_decor("K_charge");
144 SG::AuxElement::Decorator<float> ChargePi_s_decor("Pi_s_charge");
145 SG::AuxElement::Decorator<float> Chi2Mu_decor("Mu_chi2");
146 SG::AuxElement::Decorator<float> nDoFMu_decor("Mu_nDoF");
147 //muon contribution to chi2 of the cascade fit
148 SG::AuxElement::Decorator<float> MuChi2B_decor("Mu_chi2_B");
149 SG::AuxElement::Decorator<float> MunDoFB_decor("Mu_nDoF_B");
150
151
152 // Get mu+pi container and identify the input mu+pi
153 const xAOD::VertexContainer *MuPiContainer{};
154 ATH_CHECK(evtStore()->retrieve(MuPiContainer , m_vertexContainerKey )); // FIXME Use Handles
155 const xAOD::VertexContainer *d0Container{};
156 ATH_CHECK(evtStore()->retrieve(d0Container , m_vertexD0ContainerKey )); //"D0Vertices" // FIXME Use Handles
157
158 for (Trk::VxCascadeInfo* x : cascadeinfoContainer) {
159 if(x==nullptr){
160 ATH_MSG_ERROR("cascadeinfoContainer is null");
161 continue;
162 }
163
164 // the cascade fitter returns:
165 // std::vector<xAOD::Vertex*>, each xAOD::Vertex contains the refitted track parameters (perigee at the vertex position)
166 // vertices[iv] the links to the original TPs and a covariance of size 3+5*NTRK; the chi2 of the total fit
167 // is split between the cascade vertices as per track contribution
168 // std::vector< std::vector<TLorentzVector> >, each std::vector<TLorentzVector> contains the refitted momenta (TLorentzVector)
169 // momenta[iv][...] of all tracks in the corresponding vertex, including any pseudotracks (from cascade vertices)
170 // originating in this vertex; the masses are as assigned in the cascade fit
171 // std::vector<Amg::MatrixX>, the corresponding covariance matrices in momentum space
172 // covariance[iv]
173 // int nDoF, double Chi2
174 //
175 // the invariant mass, pt, lifetime etc. errors should be calculated using the covariance matrices in momentum space as these
176 // take into account the full track-track and track-vertex correlations
177 //
178 // 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.
179 // The covariance terms between the two vertices are not stored. In momentum space momenta[0] contains the 2 V0 tracks,
180 // their momenta add up to the momentum of the 3rd track in momenta[1], the first two being the Jpsi tracks
181
182 const std::vector<xAOD::Vertex*> &cascadeVertices = x->vertices();
183 if(cascadeVertices.size()!=topoN)
184 ATH_MSG_ERROR("Incorrect number of vertices");
185 if(cascadeVertices[0] == nullptr || cascadeVertices[1] == nullptr) ATH_MSG_ERROR("Error null vertex");
186 // Keep vertices (bear in mind that they come in reverse order!)
187 for(int i =0;i<topoN;i++) Vtxwritehandles[i]->push_back(cascadeVertices[i]);
188
189 x->setSVOwnership(false); // Prevent Container from deleting vertices
190 const auto mainVertex = cascadeVertices[1]; // this is the B_c+/- vertex
191 const std::vector< std::vector<TLorentzVector> > &moms = x->getParticleMoms();
192
193 // Set links to cascade vertices
194 std::vector<const xAOD::Vertex*> verticestoLink;
195 verticestoLink.push_back(cascadeVertices[0]);
196 if(Vtxwritehandles[1] == nullptr) ATH_MSG_ERROR("Vtxwritehandles[1] is null");
197 if(!BPhysPVCascadeTools::LinkVertices(CascadeLinksDecor, verticestoLink, Vtxwritehandles[0], cascadeVertices[1]))
198 ATH_MSG_ERROR("Error decorating with cascade vertices");
199
200 // Identify the input mu+pi_soft
201 const xAOD::Vertex* MuPiVertex = BPhysPVCascadeTools::FindVertex<2>(MuPiContainer, cascadeVertices[1]);
202 ATH_MSG_DEBUG("1 pt mu+pi_soft tracks " << cascadeVertices[1]->trackParticle(0)->pt() << ", " << cascadeVertices[1]->trackParticle(1)->pt());
203 if (MuPiVertex) ATH_MSG_DEBUG("2 pt mu+pi_soft tracks " << MuPiVertex->trackParticle(0)->pt() << ", " << MuPiVertex->trackParticle(1)->pt());
204
205 // Identify the input D0
206 const xAOD::Vertex* d0Vertex = BPhysPVCascadeTools::FindVertex<2>(d0Container, cascadeVertices[0]);;
207 ATH_MSG_DEBUG("1 pt D0 tracks " << cascadeVertices[0]->trackParticle(0)->pt() << ", " << cascadeVertices[0]->trackParticle(1)->pt());
208 if (d0Vertex) ATH_MSG_DEBUG("2 pt D0 tracks " << d0Vertex->trackParticle(0)->pt() << ", " << d0Vertex->trackParticle(1)->pt());
209
210 // Set links to input vertices
211 std::vector<const xAOD::Vertex*> MuPiVerticestoLink;
212 if (MuPiVertex) MuPiVerticestoLink.push_back(MuPiVertex);
213 else ATH_MSG_WARNING("Could not find linking mu+pi_soft");
214 if(!BPhysPVCascadeTools::LinkVertices(MuPiLinksDecor, MuPiVerticestoLink, MuPiContainer, cascadeVertices[1]))
215 ATH_MSG_ERROR("Error decorating with mu+pi_soft vertices");
216
217 std::vector<const xAOD::Vertex*> d0VerticestoLink;
218 if (d0Vertex) d0VerticestoLink.push_back(d0Vertex);
219 else ATH_MSG_WARNING("Could not find linking D0");
220 if(!BPhysPVCascadeTools::LinkVertices(D0LinksDecor, d0VerticestoLink, d0Container, cascadeVertices[1]))
221 ATH_MSG_ERROR("Error decorating with D0 vertices");
222
223 bool tagD0(true);
224 if (MuPiVertex){
225 if(std::abs(m_Dx_pid)==421 && (MuPiVertex->trackParticle(1)->charge()==-1)) tagD0 = false; //checking soft pion charge
226 }
227
228 double mass_b = m_vtx0MassHypo;
229 double mass_d0 = m_vtx1MassHypo;
230 std::vector<double> massesMuPi;
231 massesMuPi.push_back(m_vtx0Daug1MassHypo); //mu
232 massesMuPi.push_back(m_vtx0Daug2MassHypo); //pi_soft
233 std::vector<double> massesD0;
234 if(tagD0){
235 massesD0.push_back(m_vtx1Daug1MassHypo); //pi
236 massesD0.push_back(m_vtx1Daug2MassHypo); //K
237 }else{ // Change the order of masses for D*-->D0bar pi-, D0bar->K+pi-
238 massesD0.push_back(m_vtx1Daug2MassHypo); //K
239 massesD0.push_back(m_vtx1Daug1MassHypo); //pi
240 }
241 std::vector<double> Masses; // masses of all particles "inside" of B
242 Masses.push_back(m_vtx0Daug1MassHypo); //mu
243 Masses.push_back(m_vtx0Daug2MassHypo); //pi_soft
244 Masses.push_back(m_vtx1MassHypo); //D0 mass
245
246 // loop over candidates -- Don't apply PV_minNTracks requirement here
247 // because it may result in exclusion of the high-pt PV.
248 // get good PVs
249
250 xAOD::BPhysHypoHelper vtx(m_hypoName, mainVertex);
251
252 // Get refitted track momenta from all vertices, charged tracks only
253 BPhysPVCascadeTools::SetVectorInfo(vtx, x);
254
255
256 // Decorate main vertex
257 //
258 // 1.a) mass, mass error
259 BPHYS_CHECK( vtx.setMass(m_CascadeTools->invariantMass(moms[1])) );
260 BPHYS_CHECK( vtx.setMassErr(m_CascadeTools->invariantMassError(moms[1],x->getCovariance()[1])) );
261 // 1.b) pt and pT error (the default pt of mainVertex is != the pt of the full cascade fit!)
262 Pt_decor(*mainVertex) = m_CascadeTools->pT(moms[1]);
263 PtErr_decor(*mainVertex) = m_CascadeTools->pTError(moms[1],x->getCovariance()[1]);
264 // 1.c) chi2 and ndof (the default chi2 of mainVertex is != the chi2 of the full cascade fit!)
265 chi2_decor(*mainVertex) = x->fitChi2();
266 ndof_decor(*mainVertex) = x->nDoF();
267 Chi2Mu_decor(*mainVertex) = cascadeVertices[1]->trackParticle(0)->chiSquared();
268 nDoFMu_decor(*mainVertex) = cascadeVertices[1]->trackParticle(0)->numberDoF();
269
270 // track contribution to chi2 of cascasde fit
271 std::vector< Trk::VxTrackAtVertex > trkAtB = cascadeVertices[1]->vxTrackAtVertex();
272 MuChi2B_decor(*mainVertex) = trkAtB.at(0).trackQuality().chiSquared();
273 MunDoFB_decor(*mainVertex) = trkAtB.at(0).trackQuality().numberDoF();
274
275 float massMuPi = 0.;
276 if (MuPiVertex) {
277 TLorentzVector p4_mu1, p4_pi_s;
278 p4_mu1.SetPtEtaPhiM(MuPiVertex->trackParticle(0)->pt(),
279 MuPiVertex->trackParticle(0)->eta(),
280 MuPiVertex->trackParticle(0)->phi(), m_vtx0Daug1MassHypo); //mu
281 p4_pi_s.SetPtEtaPhiM(MuPiVertex->trackParticle(1)->pt(),
282 MuPiVertex->trackParticle(1)->eta(),
283 MuPiVertex->trackParticle(1)->phi(), m_vtx0Daug2MassHypo); // pi_soft
284 massMuPi = (p4_mu1 + p4_pi_s).M();
285 ChargePi_s_decor(*mainVertex) = MuPiVertex->trackParticle(1)->charge();
286 ChargeMu_decor(*mainVertex) = MuPiVertex->trackParticle(0)->charge();
287 }
288 massMuPi_decor(*mainVertex) = massMuPi;
289
290 //float massKpi = 0.;
291 if (d0Vertex) {
292 TLorentzVector p4_ka, p4_pi;
293 if(tagD0){
294 p4_pi.SetPtEtaPhiM(d0Vertex->trackParticle(0)->pt(),
295 d0Vertex->trackParticle(0)->eta(),
296 d0Vertex->trackParticle(0)->phi(), m_vtx1Daug1MassHypo);
297 p4_ka.SetPtEtaPhiM(d0Vertex->trackParticle(1)->pt(),
298 d0Vertex->trackParticle(1)->eta(),
299 d0Vertex->trackParticle(1)->phi(), m_vtx1Daug2MassHypo);
300 ChargeK_decor(*mainVertex) = d0Vertex->trackParticle(1)->charge();
301 ChargePi_decor(*mainVertex) = d0Vertex->trackParticle(0)->charge();
302
303 }else{ // Change the oreder of masses for D*-->D0bar pi-, D0bar->K+pi-
304 p4_pi.SetPtEtaPhiM(d0Vertex->trackParticle(1)->pt(),
305 d0Vertex->trackParticle(1)->eta(),
306 d0Vertex->trackParticle(1)->phi(), m_vtx1Daug1MassHypo);
307 p4_ka.SetPtEtaPhiM(d0Vertex->trackParticle(0)->pt(),
308 d0Vertex->trackParticle(0)->eta(),
309 d0Vertex->trackParticle(0)->phi(), m_vtx1Daug2MassHypo);
310 ChargeK_decor(*mainVertex) = d0Vertex->trackParticle(0)->charge();
311 ChargePi_decor(*mainVertex) = d0Vertex->trackParticle(1)->charge();
312 }
313 //massKpi = (p4_ka + p4_pi).M();
314 }
315 MassMuPiAft_decor(*mainVertex) = (moms[1][0] + moms[1][1]).M(); //mu + pi_soft
316 MassPiD0_decor(*mainVertex) = (moms[1][1] + moms[0][0] + moms[0][1]).M(); //pi_soft + D0
317 MassKpi_svdecor(*mainVertex) = (moms[0][0] + moms[0][1]).M();//D0
318 MassMuPiPiK_decor(*mainVertex) = m_CascadeTools->invariantMass(moms[1]); //B0
319
320 ATH_CHECK(helper.FillCandwithRefittedVertices(m_refitPV, pvContainer, refPvContainer, &(*m_pvRefitter), m_PV_max, m_DoVertexType, x, 1, mass_b, vtx));
321
322 // 4) decorate the main vertex with D0 vertex mass, pt, lifetime and lxy values (plus errors)
323 // D0 points to the main vertex, so lifetime and lxy are w.r.t the main vertex
324 Mass_svdecor(*mainVertex) = m_CascadeTools->invariantMass(moms[0]);
325 MassErr_svdecor(*mainVertex) = m_CascadeTools->invariantMassError(moms[0],x->getCovariance()[0]);
326 Pt_svdecor(*mainVertex) = m_CascadeTools->pT(moms[0]);
327 PtErr_svdecor(*mainVertex) = m_CascadeTools->pTError(moms[0],x->getCovariance()[0]);
328 Lxy_svdecor(*mainVertex) = m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]);
329 LxyErr_svdecor(*mainVertex) = m_CascadeTools->lxyError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
330 Tau_svdecor(*mainVertex) = m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1]);
331 TauErr_svdecor(*mainVertex) = m_CascadeTools->tauError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1]);
332
333
334 // Some checks in DEBUG mode
335 ATH_MSG_DEBUG("chi2 " << x->fitChi2() //DEBUG->INFO
336 << " chi2_1 " << m_V0Tools->chisq(cascadeVertices[0])
337 << " chi2_2 " << m_V0Tools->chisq(cascadeVertices[1])
338 << " vprob " << m_CascadeTools->vertexProbability(x->nDoF(),x->fitChi2()));
339 ATH_MSG_DEBUG("ndf " << x->nDoF() << " ndf_1 " << m_V0Tools->ndof(cascadeVertices[0]) << " ndf_2 " << m_V0Tools->ndof(cascadeVertices[1]));
340 ATH_MSG_DEBUG("V0Tools mass_d0 " << m_V0Tools->invariantMass(cascadeVertices[0],massesD0)
341 << " error " << m_V0Tools->invariantMassError(cascadeVertices[0],massesD0)
342 << " mass_J " << m_V0Tools->invariantMass(cascadeVertices[1],massesMuPi)
343 << " error " << m_V0Tools->invariantMassError(cascadeVertices[1],massesMuPi));
344 // masses and errors, using track masses assigned in the fit
345 double Mass_B = m_CascadeTools->invariantMass(moms[1]);
346 double Mass_D0 = m_CascadeTools->invariantMass(moms[0]);
347 double Mass_B_err = m_CascadeTools->invariantMassError(moms[1],x->getCovariance()[1]);
348 double Mass_D0_err = m_CascadeTools->invariantMassError(moms[0],x->getCovariance()[0]);
349 ATH_MSG_DEBUG("Mass_B " << Mass_B << " Mass_D0 " << Mass_D0);
350 ATH_MSG_DEBUG("Mass_B_err " << Mass_B_err << " Mass_D0_err " << Mass_D0_err);
351 double mprob_B = m_CascadeTools->massProbability(mass_b,Mass_B,Mass_B_err);
352 double mprob_D0 = m_CascadeTools->massProbability(mass_d0,Mass_D0,Mass_D0_err);
353 ATH_MSG_DEBUG("mprob_B " << mprob_B << " mprob_D0 " << mprob_D0);
354 // masses and errors, assigning user defined track masses
355 ATH_MSG_DEBUG("Mass_b " << m_CascadeTools->invariantMass(moms[1],Masses)
356 << " Mass_d0 " << m_CascadeTools->invariantMass(moms[0],massesD0));
357 ATH_MSG_DEBUG("Mass_b_err " << m_CascadeTools->invariantMassError(moms[1],x->getCovariance()[1],Masses)
358 << " Mass_d0_err " << m_CascadeTools->invariantMassError(moms[0],x->getCovariance()[0],massesD0));
359 ATH_MSG_DEBUG("pt_b " << m_CascadeTools->pT(moms[1])
360 << " pt_d " << m_CascadeTools->pT(moms[0])
361 << " pt_d0 " << m_V0Tools->pT(cascadeVertices[0]));
362 ATH_MSG_DEBUG("ptErr_b " << m_CascadeTools->pTError(moms[1],x->getCovariance()[1])
363 << " ptErr_d " << m_CascadeTools->pTError(moms[0],x->getCovariance()[0])
364 << " ptErr_d0 " << m_V0Tools->pTError(cascadeVertices[0]));
365 ATH_MSG_DEBUG("lxy_B " << m_V0Tools->lxy(cascadeVertices[1],primaryVertex) << " lxy_D " << m_V0Tools->lxy(cascadeVertices[0],cascadeVertices[1]));
366 ATH_MSG_DEBUG("lxy_b " << m_CascadeTools->lxy(moms[1],cascadeVertices[1],primaryVertex) << " lxy_d " << m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]));
367 ATH_MSG_DEBUG("lxyErr_b " << m_CascadeTools->lxyError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
368 << " lxyErr_d " << m_CascadeTools->lxyError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
369 << " lxyErr_d0 " << m_V0Tools->lxyError(cascadeVertices[0],cascadeVertices[1]));
370 ATH_MSG_DEBUG("tau_B " << m_CascadeTools->tau(moms[1],cascadeVertices[1],primaryVertex,mass_b)
371 << " tau_d0 " << m_V0Tools->tau(cascadeVertices[0],cascadeVertices[1],massesD0));
372 ATH_MSG_DEBUG("tau_b " << m_CascadeTools->tau(moms[1],cascadeVertices[1],primaryVertex)
373 << " tau_d " << m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1])
374 << " tau_D " << m_CascadeTools->tau(moms[0],cascadeVertices[0],cascadeVertices[1],mass_d0));
375 ATH_MSG_DEBUG("tauErr_b " << m_CascadeTools->tauError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
376 << " tauErr_d " << m_CascadeTools->tauError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
377 << " tauErr_d0 " << m_V0Tools->tauError(cascadeVertices[0],cascadeVertices[1],massesD0));
378 ATH_MSG_DEBUG("TauErr_b " << m_CascadeTools->tauError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex,mass_b)
379 << " TauErr_d " << m_CascadeTools->tauError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1],mass_d0)
380 << " TauErr_d0 " << m_V0Tools->tauError(cascadeVertices[0],cascadeVertices[1],massesD0,mass_d0));
381
382 ATH_MSG_DEBUG("CascadeTools main vert wrt PV " << " CascadeTools SV " << " V0Tools SV");
383 ATH_MSG_DEBUG("a0z " << m_CascadeTools->a0z(moms[1],cascadeVertices[1],primaryVertex)
384 << ", " << m_CascadeTools->a0z(moms[0],cascadeVertices[0],cascadeVertices[1])
385 << ", " << m_V0Tools->a0z(cascadeVertices[0],cascadeVertices[1]));
386 ATH_MSG_DEBUG("a0zErr " << m_CascadeTools->a0zError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
387 << ", " << m_CascadeTools->a0zError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
388 << ", " << m_V0Tools->a0zError(cascadeVertices[0],cascadeVertices[1]));
389 ATH_MSG_DEBUG("a0xy " << m_CascadeTools->a0xy(moms[1],cascadeVertices[1],primaryVertex)
390 << ", " << m_CascadeTools->a0xy(moms[0],cascadeVertices[0],cascadeVertices[1])
391 << ", " << m_V0Tools->a0xy(cascadeVertices[0],cascadeVertices[1]));
392 ATH_MSG_DEBUG("a0xyErr " << m_CascadeTools->a0xyError(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
393 << ", " << m_CascadeTools->a0xyError(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
394 << ", " << m_V0Tools->a0xyError(cascadeVertices[0],cascadeVertices[1]));
395 ATH_MSG_DEBUG("a0 " << m_CascadeTools->a0(moms[1],cascadeVertices[1],primaryVertex)
396 << ", " << m_CascadeTools->a0(moms[0],cascadeVertices[0],cascadeVertices[1])
397 << ", " << m_V0Tools->a0(cascadeVertices[0],cascadeVertices[1]));
398 ATH_MSG_DEBUG("a0Err " << m_CascadeTools->a0Error(moms[1],x->getCovariance()[1],cascadeVertices[1],primaryVertex)
399 << ", " << m_CascadeTools->a0Error(moms[0],x->getCovariance()[0],cascadeVertices[0],cascadeVertices[1])
400 << ", " << m_V0Tools->a0Error(cascadeVertices[0],cascadeVertices[1]));
401 ATH_MSG_DEBUG("x0 " << m_V0Tools->vtx(cascadeVertices[0]).x() << " y0 " << m_V0Tools->vtx(cascadeVertices[0]).y() << " z0 " << m_V0Tools->vtx(cascadeVertices[0]).z());
402 ATH_MSG_DEBUG("x1 " << m_V0Tools->vtx(cascadeVertices[1]).x() << " y1 " << m_V0Tools->vtx(cascadeVertices[1]).y() << " z1 " << m_V0Tools->vtx(cascadeVertices[1]).z());
403 ATH_MSG_DEBUG("X0 " << primaryVertex->x() << " Y0 " << primaryVertex->y() << " Z0 " << primaryVertex->z());
404 ATH_MSG_DEBUG("rxy0 " << m_V0Tools->rxy(cascadeVertices[0]) << " rxyErr0 " << m_V0Tools->rxyError(cascadeVertices[0]));
405 ATH_MSG_DEBUG("rxy1 " << m_V0Tools->rxy(cascadeVertices[1]) << " rxyErr1 " << m_V0Tools->rxyError(cascadeVertices[1]));
406 ATH_MSG_DEBUG("Rxy0 wrt PV " << m_V0Tools->rxy(cascadeVertices[0],primaryVertex) << " RxyErr0 wrt PV " << m_V0Tools->rxyError(cascadeVertices[0],primaryVertex));
407 ATH_MSG_DEBUG("Rxy1 wrt PV " << m_V0Tools->rxy(cascadeVertices[1],primaryVertex) << " RxyErr1 wrt PV " << m_V0Tools->rxyError(cascadeVertices[1],primaryVertex));
408 ATH_MSG_DEBUG("number of covariance matrices " << (x->getCovariance()).size());
409 } // loop over cascadeinfoContainer
410
411 // Deleting cascadeinfo since this won't be stored.
412 // Vertices have been kept in m_cascadeOutputs and should be owned by their container
413 for (auto x : cascadeinfoContainer) delete x;
414
415 return StatusCode::SUCCESS;
416 }
417
418
419 MuPlusDpstCascade::MuPlusDpstCascade(const std::string& t, const std::string& n, const IInterface* p) : base_class(t,n,p),
420 m_vertexContainerKey(""),
421 m_vertexD0ContainerKey(""),
422 m_cascadeOutputsKeys{ "MuPlusDpstCascadeVtx1", "MuPlusDpstCascadeVtx2" },
423 m_VxPrimaryCandidateName("PrimaryVertices"),
424 m_MuPiMassLower(0.0),
425 m_MuPiMassUpper(10000.0),
426 m_D0MassLower(0.0),
427 m_D0MassUpper(10000.0),
428 m_DstMassLower(0.0),
429 m_DstMassUpper(10000.0),
430 m_DstMassUpperAft(10000.0),
431 m_MassLower(0.0),
432 m_MassUpper(20000.0),
433 m_vtx0MassHypo(-1),
434 m_vtx1MassHypo(-1),
435 m_vtx0Daug1MassHypo(-1),
436 m_vtx0Daug2MassHypo(-1),
437 m_vtx1Daug1MassHypo(-1),
438 m_vtx1Daug2MassHypo(-1),
439 m_Dx_pid(421),
440 m_constrD0(true),
441 m_constrMuPi(false),
442 m_chi2cut(-1.0),
443 m_iVertexFitter("Trk::TrkVKalVrtFitter"),
444 m_pvRefitter("Analysis::PrimaryVertexRefitter"),
445 m_V0Tools("Trk::V0Tools"),
446 m_CascadeTools("DerivationFramework::CascadeTools")
447 {
448 declareProperty("MuPiVertices", m_vertexContainerKey);
449 declareProperty("D0Vertices", m_vertexD0ContainerKey);
450 declareProperty("VxPrimaryCandidateName", m_VxPrimaryCandidateName);
451 declareProperty("RefPVContainerName", m_refPVContainerName = "RefittedPrimaryVertices");
452 declareProperty("MuPiMassLowerCut", m_MuPiMassLower);
453 declareProperty("MuPiMassUpperCut", m_MuPiMassUpper);
454 declareProperty("D0MassLowerCut", m_D0MassLower);
455 declareProperty("D0MassUpperCut", m_D0MassUpper);
456 declareProperty("DstMassLowerCut", m_DstMassLower);
457 declareProperty("DstMassUpperCut", m_DstMassUpper);
458 declareProperty("DstMassUpperCutAft", m_DstMassUpperAft); //mass cut after cascade fit
459 declareProperty("MassLowerCut", m_MassLower);
460 declareProperty("MassUpperCut", m_MassUpper);
461 declareProperty("HypothesisName", m_hypoName = "B");
462 declareProperty("Vtx0MassHypo", m_vtx0MassHypo);
463 declareProperty("Vtx1MassHypo", m_vtx1MassHypo);
464 declareProperty("Vtx0Daug1MassHypo", m_vtx0Daug1MassHypo);
465 declareProperty("Vtx0Daug2MassHypo", m_vtx0Daug2MassHypo);
466 declareProperty("Vtx0Daug3MassHypo", m_vtx0Daug2MassHypo);
467 declareProperty("Vtx1Daug1MassHypo", m_vtx1Daug1MassHypo);
468 declareProperty("Vtx1Daug2MassHypo", m_vtx1Daug2MassHypo);
469 declareProperty("DxHypothesis", m_Dx_pid);
470 declareProperty("ApplyD0MassConstraint", m_constrD0);
471 declareProperty("ApplyMuPiMassConstraint", m_constrMuPi);
472 declareProperty("Chi2Cut", m_chi2cut);
473 declareProperty("RefitPV", m_refitPV = true);
474 declareProperty("MaxnPV", m_PV_max = 999);
475 declareProperty("MinNTracksInPV", m_PV_minNTracks = 0);
476 declareProperty("DoVertexType", m_DoVertexType = 7);
477 declareProperty("TrkVertexFitterTool", m_iVertexFitter);
478 declareProperty("PVRefitter", m_pvRefitter);
479 declareProperty("V0Tools", m_V0Tools);
480 declareProperty("CascadeTools", m_CascadeTools);
481 declareProperty("CascadeVertexCollections", m_cascadeOutputsKeys);
482 }
483
484 MuPlusDpstCascade::~MuPlusDpstCascade(){ }
485
486 StatusCode MuPlusDpstCascade::performSearch(std::vector<Trk::VxCascadeInfo*> *cascadeinfoContainer, const EventContext& ctx) const
487 {
488
489 assert(cascadeinfoContainer!=nullptr);
490
491 // Get TrackParticle container (for setting links to the original tracks)
492 const xAOD::TrackParticleContainer *trackContainer{};
493 ATH_CHECK(evtStore()->retrieve(trackContainer , "InDetTrackParticles" )); // FIXME Use Handles
494
495 // Get mu+pi container
496 const xAOD::VertexContainer *MuPiContainer{};
497 ATH_CHECK(evtStore()->retrieve(MuPiContainer , m_vertexContainerKey )); // FIXME Use Handles
498
499 // Get D0 container
500 const xAOD::VertexContainer *d0Container{};
501 ATH_CHECK(evtStore()->retrieve(d0Container , m_vertexD0ContainerKey )); //"D0Vertices" // FIXME Use Handles
502
503 double mass_d0 = m_vtx1MassHypo;
504 std::vector<const xAOD::TrackParticle*> tracksMuPi;
505 std::vector<const xAOD::TrackParticle*> tracksD0;
506 std::vector<double> massesMuPi;
507 massesMuPi.push_back(m_vtx0Daug1MassHypo); //mass mu
508 massesMuPi.push_back(m_vtx0Daug2MassHypo); //mass pi
509 std::vector<double> massesD0;
510 massesD0.push_back(m_vtx1Daug1MassHypo); //mass pi
511 massesD0.push_back(m_vtx1Daug2MassHypo); //mass K
512 std::vector<double> massesD0b; // Change the oreder of masses for D*-->D0bar pi-, D0bar->K+pi-
513 massesD0b.push_back(m_vtx1Daug2MassHypo);
514 massesD0b.push_back(m_vtx1Daug1MassHypo);
515 std::vector<double> Masses;
516 Masses.push_back(m_vtx0Daug1MassHypo); //mu
517 Masses.push_back(m_vtx0Daug2MassHypo); //pi
518 Masses.push_back(m_vtx1MassHypo); //D0
519
520
521
522 // Select mu+pi_soft candidates before calling cascade fit
523 std::vector<const xAOD::Vertex*> selectedMuPiCandidates;
524 for ( auto vxcItr : *MuPiContainer ){
525 // Check mu+pi_soft candidate invariant mass and skip if need be
526 TLorentzVector p4Mup_in, p4Mum_in;
527 p4Mup_in.SetPtEtaPhiM(vxcItr->trackParticle(0)->pt(),
528 vxcItr->trackParticle(0)->eta(),
529 vxcItr->trackParticle(0)->phi(), m_vtx0Daug1MassHypo);
530 p4Mum_in.SetPtEtaPhiM(vxcItr->trackParticle(1)->pt(),
531 vxcItr->trackParticle(1)->eta(),
532 vxcItr->trackParticle(1)->phi(), m_vtx0Daug2MassHypo);
533 double mass_MuPi = (p4Mup_in + p4Mum_in).M();
534 ATH_MSG_DEBUG("mu pi_soft mass " << mass_MuPi);
535 if (mass_MuPi < m_MuPiMassLower || mass_MuPi > m_MuPiMassUpper) {
536 ATH_MSG_DEBUG(" Original mu & pi_soft candidate rejected by the mass cut: mass = "
537 << mass_MuPi << " != (" << m_MuPiMassLower << ", " << m_MuPiMassUpper << ")" );
538 continue;
539 }
540
541 // Track selection - Loose
542 // for soft pion wich is (2nd) in MuPi vertex
543 if ( !m_trackSelectionTools->accept(vxcItr->trackParticle(1)) ){
544 ATH_MSG_DEBUG(" Original mu & pi_soft candidate rejected by the track's cut level - loose");
545 continue;
546 }
547
548 selectedMuPiCandidates.push_back(vxcItr);
549 } //for(auto vxcItr : *MuPiContainer)
550 if(selectedMuPiCandidates.size()<1) return StatusCode::SUCCESS;
551
552 // Select the D0/D0b candidates before calling cascade fit
553 std::vector<const xAOD::Vertex*> selectedD0Candidates;
554 for(auto vxcItr : *d0Container){
555 // Check the passed flag first
556 const xAOD::Vertex* vtx = vxcItr;
557 SG::AuxElement::Accessor<Char_t> flagAcc1("passed_D0");
558 SG::AuxElement::Accessor<Char_t> flagAcc2("passed_D0b");
559 bool isD0(true);
560 bool isD0b(true);
561 if(flagAcc1.isAvailable(*vtx)){
562 if(!flagAcc1(*vtx)) isD0 = false;
563 }
564 if(flagAcc2.isAvailable(*vtx)){
565 if(!flagAcc2(*vtx)) isD0b = false;
566 }
567 if(!(isD0||isD0b)) continue;
568
569 // Track selection - Loose
570 if ( !m_trackSelectionTools->accept(vxcItr->trackParticle(0)) ){
571 ATH_MSG_DEBUG(" Original D0/D0-bar candidate rejected by the track's cut level - loose ");
572 continue;
573 }
574 if ( !m_trackSelectionTools->accept(vxcItr->trackParticle(1)) ){
575 ATH_MSG_DEBUG(" Original D0/D0-bar candidate rejected by the track's cut level - loose ");
576 continue;
577 }
578
579
580 // Ensure the total charge is correct
581 if (vxcItr->trackParticle(0)->charge() != 1 || vxcItr->trackParticle(1)->charge() != -1) {
582 ATH_MSG_DEBUG(" Original D0/D0-bar candidate rejected by the charge requirement: "
583 << vxcItr->trackParticle(0)->charge() << ", " << vxcItr->trackParticle(1)->charge() );
584 continue;
585 }
586
587 // Check D0/D0bar candidate invariant mass and skip if need be
588 double mass_D0 = m_V0Tools->invariantMass(vxcItr,massesD0);
589 double mass_D0b = m_V0Tools->invariantMass(vxcItr,massesD0b);
590 ATH_MSG_DEBUG("D0 mass " << mass_D0 << ", D0b mass "<<mass_D0b);
591 if ((mass_D0 < m_D0MassLower || mass_D0 > m_D0MassUpper) && (mass_D0b < m_D0MassLower || mass_D0b > m_D0MassUpper)) {
592 ATH_MSG_DEBUG(" Original D0/D0-bar candidate rejected by the mass cut: mass = "
593 << mass_D0 << " != (" << m_D0MassLower << ", " << m_D0MassUpper << ") "
594 << mass_D0b << " != (" << m_D0MassLower << ", " << m_D0MassUpper << ") " );
595 continue;
596 }
597
598 selectedD0Candidates.push_back(vxcItr);
599 } //for(auto vxcItr : *d0Container)
600 if(selectedD0Candidates.size()<1) return StatusCode::SUCCESS;
601
602 // Select mu D*+ candidates
603 // Iterate over mu+pi_soft vertices
604 for(auto MuPiItr:selectedMuPiCandidates){
605 size_t MuPiTrkNum = MuPiItr->nTrackParticles();
606
607 tracksMuPi.clear();
608 for( unsigned int it=0; it<MuPiTrkNum; it++) tracksMuPi.push_back(MuPiItr->trackParticle(it));
609
610 if (tracksMuPi.size() != 2 || massesMuPi.size() != 2 ) {
611 ATH_MSG_INFO("problems with mu+pi_soft input");
612 }
613
614 bool tagD0(true);
615 if(std::abs(m_Dx_pid)==421 && MuPiItr->trackParticle(1)->charge()==-1) tagD0 = false;
616
617 TLorentzVector p4_pi1; // Momentum of soft pion1 = our soft pion (1)
618 p4_pi1.SetPtEtaPhiM(MuPiItr->trackParticle(1)->pt(),
619 MuPiItr->trackParticle(1)->eta(),
620 MuPiItr->trackParticle(1)->phi(), m_vtx0Daug2MassHypo);
621 // Iterate over D0/D0bar vertices
622 for(auto d0Itr : selectedD0Candidates){
623
624 // Check identical tracks in input
625 if(std::find(tracksMuPi.cbegin(), tracksMuPi.cend(), d0Itr->trackParticle(0)) != tracksMuPi.cend()) continue;
626 if(std::find(tracksMuPi.cbegin(), tracksMuPi.cend(), d0Itr->trackParticle(1)) != tracksMuPi.cend()) continue;
627
628 TLorentzVector p4_ka, p4_pi2;
629 if(tagD0){ // for D*+
630 p4_pi2.SetPtEtaPhiM(d0Itr->trackParticle(0)->pt(),
631 d0Itr->trackParticle(0)->eta(),
632 d0Itr->trackParticle(0)->phi(), m_vtx1Daug1MassHypo);
633 p4_ka.SetPtEtaPhiM( d0Itr->trackParticle(1)->pt(),
634 d0Itr->trackParticle(1)->eta(),
635 d0Itr->trackParticle(1)->phi(), m_vtx1Daug2MassHypo);
636 }else{ // change the order in the case of D*-
637 p4_pi2.SetPtEtaPhiM(d0Itr->trackParticle(1)->pt(),
638 d0Itr->trackParticle(1)->eta(),
639 d0Itr->trackParticle(1)->phi(), m_vtx1Daug1MassHypo);
640 p4_ka.SetPtEtaPhiM( d0Itr->trackParticle(0)->pt(),
641 d0Itr->trackParticle(0)->eta(),
642 d0Itr->trackParticle(0)->phi(), m_vtx1Daug2MassHypo);
643 }
644 // Check D*+/- candidate invariant mass and skip if need be
645 double mass_Dst= (p4_pi1 + p4_ka + p4_pi2).M();
646 ATH_MSG_DEBUG("D*+/- mass " << mass_Dst);
647 if (mass_Dst < m_DstMassLower || mass_Dst > m_DstMassUpper) {
648 ATH_MSG_DEBUG(" Original D*+/- candidate rejected by the mass cut: mass = "
649 << mass_Dst << " != (" << m_DstMassLower << ", " << m_DstMassUpper << ")" );
650 continue;
651 }
652
653 size_t d0TrkNum = d0Itr->nTrackParticles(); //2
654 tracksD0.clear();
655 for( unsigned int it=0; it<d0TrkNum; it++) tracksD0.push_back(d0Itr->trackParticle(it));
656 if (tracksD0.size() != 2 || massesD0.size() != 2 ) {
657 ATH_MSG_INFO("problems with D0 input");
658 }
659
660 //Leaving for the possible mods in the future
661 //SG::AuxElement::Accessor<xAOD::Vertex_v1::TrackParticleLinks_t> trackAcc( "trackParticleLinks" );
662 //ATH_MSG_INFO("CUSTOM:: "<<*(trackAcc(*d0Itr)).at(0));
663 //ATH_MSG_INFO("CUSTOM2:: "<<tracksD0.at(0)); //-> gives the same result
664
665 // Apply the user's settings to the fitter
666 // Reset
667 std::unique_ptr<Trk::IVKalState> state (m_iVertexFitter->makeState(ctx));
668 // Robustness
669 int robustness = 0;
670 m_iVertexFitter->setRobustness(robustness, *state);
671 // Build up the topology
672 // Vertex list
673 std::vector<Trk::VertexID> vrtList;
674 // D0 vertex
675 Trk::VertexID vID;
676 if (m_constrD0) { //ApplyD0MassConstraint = true
677 if(tagD0) vID = m_iVertexFitter->startVertex(tracksD0,massesD0, *state,mass_d0);
678 else vID = m_iVertexFitter->startVertex(tracksD0,massesD0b, *state,mass_d0);
679 } else {
680 if(tagD0) vID = m_iVertexFitter->startVertex(tracksD0, massesD0,*state);
681 else vID = m_iVertexFitter->startVertex(tracksD0, massesD0b, *state);
682 }
683 vrtList.push_back(vID);
684 // B vertex including mu+pi_soft
685 Trk::VertexID vID2 = m_iVertexFitter->nextVertex(tracksMuPi,massesMuPi,vrtList,*state);
686 if (m_constrMuPi) {
687 std::vector<Trk::VertexID> cnstV;
688 cnstV.clear();
689 if ( !m_iVertexFitter->addMassConstraint(vID2,tracksMuPi,cnstV,*state,m_vtx0MassHypo).isSuccess() ) {
690 ATH_MSG_WARNING("addMassConstraint failed");
691 //return StatusCode::FAILURE;
692 }
693 }
694
695 // Do the work
696 std::unique_ptr<Trk::VxCascadeInfo> result(m_iVertexFitter->fitCascade(*state));
697
698 if (result != nullptr) {
699 // reset links to original tracks
700 BPhysPVCascadeTools::PrepareVertexLinks(result.get(), trackContainer);
701 ATH_MSG_DEBUG("storing tracks " << ((result->vertices())[0])->trackParticle(0) << ", "
702 << ((result->vertices())[0])->trackParticle(1) << ", "
703 << ((result->vertices())[1])->trackParticle(0) << ", "
704 << ((result->vertices())[1])->trackParticle(1));
705 // necessary to prevent memory leak
706 result->setSVOwnership(true);
707
708 // Chi2/DOF cut
709 double bChi2DOF = result->fitChi2()/result->nDoF();
710 ATH_MSG_DEBUG("Candidate chi2/DOF is " << bChi2DOF);
711 bool chi2CutPassed = (m_chi2cut <= 0.0 || bChi2DOF < m_chi2cut);
712
713 const std::vector< std::vector<TLorentzVector> > &moms = result->getParticleMoms();
714 const std::vector<xAOD::Vertex*> &cascadeVertices = result->vertices();
715
716 double mass = m_CascadeTools->invariantMass(moms[1]);
717 double DstMassAft = (moms[1][1] + moms[0][0] + moms[0][1]).M(); //pi_soft + D0
718
719 if(chi2CutPassed) {
720 if (mass >= m_MassLower && mass <= m_MassUpper) {
721 if (m_CascadeTools->pT(moms[1]) > 9500){ //B_pT
722 if (m_CascadeTools->lxy(moms[0],cascadeVertices[0],cascadeVertices[1]) > 0){ //D0_Lxy>0
723 if (DstMassAft < m_DstMassUpperAft){
724
725 cascadeinfoContainer->push_back(result.release());
726
727 } //Dst_m < m_DstMassUpperAft
728 } //D0_Lxy>0
729 } //B_pT
730 } else {
731 ATH_MSG_DEBUG("Candidate rejected by the mass cut: mass = "
732 << mass << " != (" << m_MassLower << ", " << m_MassUpper << ")" );
733 } //mass Upper/Lower
734 } //chi2CutPassed
735 } //if (result != nullptr)
736
737 } //Iterate over D0 vertices
738
739 } //Iterate over mu+pi_soft vertices
740 ATH_MSG_DEBUG("cascadeinfoContainer size " << cascadeinfoContainer->size());
741 return StatusCode::SUCCESS;
742 }
743
744}
745
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.
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.
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::MuPlusDpstCascade::m_VxPrimaryCandidateName
std::string m_VxPrimaryCandidateName
Name of primary vertex container // FIXME Use Handles.
Definition MuPlusDpstCascade.h:51
DerivationFramework::MuPlusDpstCascade::m_CascadeTools
PublicToolHandle< DerivationFramework::CascadeTools > m_CascadeTools
Definition MuPlusDpstCascade.h:79
DerivationFramework::MuPlusDpstCascade::initialize
virtual StatusCode initialize() override
Definition MuPlusDpstCascade.cxx:28
DerivationFramework::MuPlusDpstCascade::m_pvRefitter
PublicToolHandle< Analysis::PrimaryVertexRefitter > m_pvRefitter
Definition MuPlusDpstCascade.h:77
DerivationFramework::MuPlusDpstCascade::m_eventInfo_key
SG::ReadHandleKey< xAOD::EventInfo > m_eventInfo_key
Definition MuPlusDpstCascade.h:75
DerivationFramework::MuPlusDpstCascade::MuPlusDpstCascade
MuPlusDpstCascade(const std::string &t, const std::string &n, const IInterface *p)
Definition MuPlusDpstCascade.cxx:419
DerivationFramework::MuPlusDpstCascade::m_V0Tools
PublicToolHandle< Trk::V0Tools > m_V0Tools
Definition MuPlusDpstCascade.h:78
DerivationFramework::MuPlusDpstCascade::m_cascadeOutputsKeys
std::vector< std::string > m_cascadeOutputsKeys
Definition MuPlusDpstCascade.h:49
DerivationFramework::MuPlusDpstCascade::m_hypoName
std::string m_hypoName
name of the mass hypothesis.
Definition MuPlusDpstCascade.h:84
DerivationFramework::MuPlusDpstCascade::performSearch
StatusCode performSearch(std::vector< Trk::VxCascadeInfo * > *cascadeinfoContainer, const EventContext &ctx) const
Definition MuPlusDpstCascade.cxx:486
DerivationFramework::MuPlusDpstCascade::m_trackSelectionTools
std::unique_ptr< InDet::InDetTrackSelectionTool > m_trackSelectionTools
Definition MuPlusDpstCascade.h:80
DerivationFramework::MuPlusDpstCascade::addBranches
virtual StatusCode addBranches(const EventContext &ctx) const override
Definition MuPlusDpstCascade.cxx:58
DerivationFramework::MuPlusDpstCascade::m_iVertexFitter
ToolHandle< Trk::TrkVKalVrtFitter > m_iVertexFitter
Definition MuPlusDpstCascade.h:76
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::TrackBag
std::vector< const xAOD::TrackParticle * > TrackBag
Definition JpsiPlusDpstCascade.cxx:25
DerivationFramework::VertexLink
ElementLink< xAOD::VertexContainer > VertexLink
Definition JpsiPlusDpstCascade.cxx:23
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::BCPLUS
static const int BCPLUS
Definition HepMCHelpers.h:124
MC::PIPLUS
static const int PIPLUS
Definition HepMCHelpers.h:110
MC::D0
static const int D0
Definition HepMCHelpers.h:120
xAOD::VertexAuxContainer
VertexAuxContainer_v1 VertexAuxContainer
Definition of the current jet auxiliary container.
Definition VertexAuxContainer.h:19
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
Generated on for ATLAS Offline Software by doxygen 1.16.1