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