ATLAS Offline Software
BPhysConversionFinder.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
5 // BPhysConversionFinder.cxx, (c) ATLAS Detector software
7 // Author: A. Chisholm <andrew.chisholm@cern.ch>
8 #include "DerivationFrameworkBPhys/BPhysConversionFinder.h"
9 #include "xAODTracking/VertexContainer.h"
10 #include "xAODTracking/VertexAuxContainer.h"
11 #include "TrkVertexAnalysisUtils/V0Tools.h"
12 #include "GaudiKernel/IPartPropSvc.h"
13 #include "GeoPrimitives/GeoPrimitivesHelpers.h"
14 #include "TrkVertexFitterInterfaces/IVertexFitter.h"
15 #include "TrkVKalVrtFitter/TrkVKalVrtFitter.h"
16 #include "TrkVKalVrtFitter/VxCascadeInfo.h"
17 #include "TLorentzVector.h"
18 #include "AthContainers/ConstAccessor.h"
19 
20 namespace DerivationFramework {
21 
22  BPhysConversionFinder::BPhysConversionFinder(const std::string& t,
23  const std::string& n,
24  const IInterface* p) :
25  AthAlgTool(t,n,p),
26  m_v0Tools("Trk::V0Tools"),
27  m_vertexFitter("Trk::TrkVKalVrtFitter"),
28  m_vertexEstimator("InDet::VertexPointEstimator"),
29  m_distanceTool("Trk::SeedNewtonTrkDistanceFinder/InDetConversionTrkDistanceFinder"),
30  m_postSelector("InDet::ConversionPostSelector"),
31  m_cascadeFitter("Trk::TrkVKalVrtFitter"),
32  m_inputTrackParticleContainerName("InDetTrackParticles"),
33  m_conversionContainerName("BPhysConversionCandidates"),
34  m_maxDistBetweenTracks(10.0),
35  m_maxDeltaCotTheta(0.3),
36  m_requireDeltaM(true),
37  m_maxDeltaM(3000.0)
38  {
39  declareInterface<DerivationFramework::IAugmentationTool>(this);
40 
41  // Declare user-defined properties
42  declareProperty("DiMuonVertexContainer", m_diMuonCollectionToCheck);
43  declareProperty("PassFlagsToCheck", m_passFlagsToCheck);
44  declareProperty("V0Tools", m_v0Tools);
45  declareProperty("VertexFitterTool", m_vertexFitter);
46  declareProperty("VertexEstimator", m_vertexEstimator);
47  declareProperty("DistanceTool", m_distanceTool);
48  declareProperty("ConversionPostSelector", m_postSelector);
49  declareProperty("CascadeFitter", m_cascadeFitter);
50  declareProperty("InputTrackParticleContainerName", m_inputTrackParticleContainerName);
51  declareProperty("ConversionContainerName", m_conversionContainerName);
52  declareProperty("MaxDistBetweenTracks", m_maxDistBetweenTracks = 10.0); // Maximum allowed distance of minimum approach
53  declareProperty("MaxDeltaCotTheta", m_maxDeltaCotTheta = 0.3); // Maximum allowed dCotTheta between tracks
54  declareProperty("RequireDeltaM", m_requireDeltaM = true); // Only save a conversions if it's a chi_c,b candidate (must then pass "MaxDeltaM" requirement), if "False" all conversions in the event will be saved
55  declareProperty("MaxDeltaM", m_maxDeltaM = 3000.0); // Maximum mass difference between di-muon+conversion and di-muon
56 
57  }
58 
59  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
60 
61  StatusCode BPhysConversionFinder::initialize()
62  {
63 
64  ATH_MSG_DEBUG("in initialize()");
65 
66  ATH_CHECK( m_v0Tools.retrieve() );
67  ATH_CHECK( m_vertexFitter.retrieve() );
68  ATH_CHECK( m_vertexEstimator.retrieve() );
69  ATH_CHECK( m_distanceTool.retrieve() );
70  ATH_CHECK( m_postSelector.retrieve() );
71  ATH_CHECK( m_cascadeFitter.retrieve() );
72 
73  return StatusCode::SUCCESS;
74 
75  }
76 
77  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
78 
79  StatusCode BPhysConversionFinder::finalize()
80  {
81  // everything all right
82  return StatusCode::SUCCESS;
83  }
84 
85  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
86 
87  StatusCode BPhysConversionFinder::addBranches() const
88  {
89 
90  int nTrackPairs_Init = 0;
91  int nTrackPairs_Selected = 0;
92  int nConv_VertexFit = 0;
93  int nConv_Selected = 0;
94  int nConv_Selected_DeltaM = 0;
95 
96  std::vector<const xAOD::Vertex*> oniaVertices;
97  oniaVertices.clear();
98 
99  //------------------------------------
100  // Look for di-muons
101  //------------------------------------
102  const xAOD::VertexContainer* diMuonContainer = NULL;
103  ATH_CHECK( evtStore()->retrieve(diMuonContainer, m_diMuonCollectionToCheck) );
104 
105  if(diMuonContainer->size() == 0) {
106 
107  ATH_MSG_DEBUG("Vertex Container (" << m_diMuonCollectionToCheck << ") is empty");
108 
109  } else {
110 
111  ATH_MSG_DEBUG("Vertex Container (" << m_diMuonCollectionToCheck << ") contains " << diMuonContainer->size() << " vertices");
112 
113  for(xAOD::VertexContainer::const_iterator vtxItr = diMuonContainer->begin(); vtxItr != diMuonContainer->end(); ++vtxItr) {
114 
115  const xAOD::Vertex* vertex = (*vtxItr);
116 
117  bool passedHypothesis = false;
118 
119  for(const auto &flag : m_passFlagsToCheck) {
120  SG::ConstAccessor<Char_t> acc(flag);
121  bool pass = acc(*vertex);
122  if(pass) passedHypothesis = true;
123  }
124 
125  if(passedHypothesis) {
126  oniaVertices.push_back(vertex);
127  }
128 
129  }
130  }
131  //------------------------------------
132 
133  // Output conversion container
134  std::unique_ptr<xAOD::VertexContainer> conversionContainer( new xAOD::VertexContainer() );
135  std::unique_ptr<xAOD::VertexAuxContainer> conversionAuxContainer( new xAOD::VertexAuxContainer() );
136  conversionContainer->setStore(conversionAuxContainer.get());
137 
138  // Only call conversion finder if we've found a di-muon candidate or
139  // we really want to look for conversions independently
140  const bool callConvFinder = !m_requireDeltaM || oniaVertices.size() > 0;
141 
142  if(callConvFinder) {
143 
144  // Retrieve track particles from StoreGate
145  const xAOD::TrackParticleContainer* inputTrackParticles = NULL;
146  ATH_CHECK( evtStore()->retrieve(inputTrackParticles,m_inputTrackParticleContainerName) );
147 
148  ATH_MSG_DEBUG("Track particle container size " << inputTrackParticles->size());
149 
150  // Track Selection
151  std::vector<const xAOD::TrackParticle*> posTracks; posTracks.clear();
152  std::vector<const xAOD::TrackParticle*> negTracks; negTracks.clear();
153 
154  // Track Loop
155  for(xAOD::TrackParticleContainer::const_iterator trkItr = inputTrackParticles->begin(); trkItr != inputTrackParticles->end(); ++trkItr) {
156 
157  const xAOD::TrackParticle* track = (*trkItr);
158 
159  uint8_t nSCT(0);
160  uint8_t nPIX(0);
161 
162  track->summaryValue(nPIX,xAOD::numberOfPixelHits);
163  track->summaryValue(nSCT,xAOD::numberOfSCTHits);
164 
165  // Don't want TRT-only tracks
166  // Require Si hits on all tracks
167  if( nSCT + nPIX < 1 ) continue;
168 
169  if( track->charge() > 0.0) {
170  posTracks.push_back(track);
171  } else {
172  negTracks.push_back(track);
173  }
174 
175  } // Track Loop
176 
177  ATH_MSG_DEBUG(posTracks.size() + negTracks.size() << " tracks pass pre-selection");
178 
179  std::vector<const xAOD::TrackParticle*>::const_iterator tpIt1;
180  std::vector<const xAOD::TrackParticle*>::const_iterator tpIt2;
181 
182  // Pos Track Loop
183  for(tpIt1 = posTracks.begin(); tpIt1 != posTracks.end(); ++tpIt1) {
184 
185  const xAOD::TrackParticle* trackParticle1 = (*tpIt1);
186 
187  const Trk::Perigee& trackPerigee1 = trackParticle1->perigeeParameters();
188 
189  // Neg Track Loop
190  for(tpIt2 = negTracks.begin(); tpIt2 != negTracks.end(); ++tpIt2) {
191 
192  if (*tpIt1 == *tpIt2) continue;
193 
194  const xAOD::TrackParticle* trackParticle2 = (*tpIt2);
195 
196  const Trk::Perigee& trackPerigee2 = trackParticle2->perigeeParameters();
197 
198  nTrackPairs_Init++;
199 
200  //------------------------------------
201  // Track pair selection
202  //------------------------------------
203  const double deltaCotTheta = fabs(1./tan(trackPerigee1.parameters()[Trk::theta]) - 1./tan(trackPerigee2.parameters()[Trk::theta]));
204  if(deltaCotTheta > m_maxDeltaCotTheta) continue;
205 
206  double distance = 1000000.;
207  std::optional<std::pair<Amg::Vector3D,Amg::Vector3D>> result = m_distanceTool->CalculateMinimumDistance(trackParticle1->perigeeParameters(),trackParticle2->perigeeParameters() );
208  bool gotDistance = result.has_value();
209  if(gotDistance) distance = Amg::distance (result->first, result->second);
210  if(!gotDistance || (distance > m_maxDistBetweenTracks)) continue;
211  //------------------------------------
212 
213  //------------------------------------
214  // Estimate starting point + cuts on compatiblity of tracks
215  //------------------------------------
216  int sflag = 0;
217  int errorcode = 0;
218  std::map<std::string, float> vertexOutput;
219  Amg::Vector3D startingPoint = m_vertexEstimator->getCirclesIntersectionPoint(&trackPerigee1,&trackPerigee2,sflag,errorcode, vertexOutput);
220  if(errorcode != 0) continue;
221  //------------------------------------
222 
223  nTrackPairs_Selected++;
224 
225  std::vector<const xAOD::TrackParticle*> trackPair;
226  trackPair.clear();
227  trackPair.push_back(trackParticle1);
228  trackPair.push_back(trackParticle2);
229 
230  // Do the vertex fit
231  std::unique_ptr<xAOD::Vertex> convVertexCandidate( m_vertexFitter->fit(trackPair, startingPoint) );
232 
233  // Check for successful fit
234  if(convVertexCandidate != NULL) {
235 
236  ATH_MSG_DEBUG("Vertex Fit Succeeded");
237 
238  convVertexCandidate->clearTracks();
239  ElementLink<xAOD::TrackParticleContainer> newLink1;
240  newLink1.setElement(*tpIt1);
241  newLink1.setStorableObject(*inputTrackParticles);
242  ElementLink<xAOD::TrackParticleContainer> newLink2;
243  newLink2.setElement(*tpIt2);
244  newLink2.setStorableObject(*inputTrackParticles);
245  convVertexCandidate->addTrackAtVertex(newLink1);
246  convVertexCandidate->addTrackAtVertex(newLink2);
247 
248  nConv_VertexFit++;
249 
250  //------------------------------------
251  // Post-vertexing cuts
252  //------------------------------------
253 
254  // This is empty and only present for compatiblity.
255  // The cut this informtion pertains to is not used for Si-Si conversions so this is OK
256  std::vector<Amg::Vector3D> positionList;
257 
258  // Apply Si-Si converion post-selection
259  if( !m_postSelector->selectConversionCandidate(convVertexCandidate.get(),0,positionList) ) {
260  convVertexCandidate.reset();
261  continue;
262  }
263  //------------------------------------
264 
265  nConv_Selected++;
266 
267  // Get photon momentum 3-vector
268  const xAOD::Vertex * constConvVertex = convVertexCandidate.get();
269  Amg::Vector3D momentum = m_v0Tools->V0Momentum(constConvVertex);
270 
271  TLorentzVector photon;
272  photon.SetXYZM(momentum.x(),momentum.y(),momentum.z(),0.0);
273 
274  //------------------------------------
275  // Check if conversion is consistent with a chi_c,b candidate
276  // by requiring a small mass difference w.r.t. any di-muon in event
277  //------------------------------------
278  bool passDeltaM = false;
279 
280  // Use to keep track of which dimuon(s) gave a chi_c/b candidate
281  std::vector<const xAOD::Vertex*> candidateOniaVertices;
282  candidateOniaVertices.clear();
283 
284  for ( std::vector<const xAOD::Vertex*>::const_iterator vtxItr = oniaVertices.begin(); vtxItr != oniaVertices.end(); ++vtxItr ) {
285 
286  const xAOD::Vertex* oniaVertex = (*vtxItr);
287 
288  static const SG::ConstAccessor<std::vector<float> > RefTrackPxAcc("RefTrackPx");
289  static const SG::ConstAccessor<std::vector<float> > RefTrackPyAcc("RefTrackPy");
290  static const SG::ConstAccessor<std::vector<float> > RefTrackPzAcc("RefTrackPz");
291  std::vector<float> diMuon_Px = RefTrackPxAcc(*oniaVertex);
292  std::vector<float> diMuon_Py = RefTrackPyAcc(*oniaVertex);
293  std::vector<float> diMuon_Pz = RefTrackPzAcc(*oniaVertex);
294 
295  TLorentzVector muon1, muon2;
296  muon1.SetXYZM(diMuon_Px.at(0),diMuon_Py.at(0),diMuon_Pz.at(0),105.658);
297  muon2.SetXYZM(diMuon_Px.at(1),diMuon_Py.at(1),diMuon_Pz.at(1),105.658);
298 
299  TLorentzVector diMuon = muon1 + muon2;
300 
301  const double deltaM = (diMuon+photon).M() - diMuon.M();
302 
303  ATH_MSG_DEBUG("Candidate DeltaM = " << deltaM << " MeV DiMuon " << oniaVertex->index() << " ( Mass = " << diMuon.M() << " MeV )");
304 
305  // Did we find a one di-muon + photon candidate with a mass diff. consistent with chi_c/b?
306  if(deltaM < m_maxDeltaM) {
307  passDeltaM = true;
308  candidateOniaVertices.push_back(oniaVertex);
309  }
310 
311  }
312 
313  // Only keep the conversion candidate if it's consistent with a chi_c,b decay
314  if(m_requireDeltaM && !passDeltaM) {
315  convVertexCandidate.reset();
316  continue;
317  }
318  //------------------------------------
319 
320  //------------------------------------
321  // Final conversion candidates
322  //------------------------------------
323  nConv_Selected_DeltaM++;
324 
325  // Keep track of which dimuon(s) gave a chi_c/b candidate
326  std::vector< ElementLink<xAOD::VertexContainer> > diMuonLinks;
327  diMuonLinks.clear();
328 
329  // Output of cascade fits with various di-muon mass hypotheses
330  std::vector<float> fit_Psi1S_Px, fit_Psi1S_Py, fit_Psi1S_Pz, fit_Psi1S_M, fit_Psi1S_ChiSq;
331  std::vector<float> fit_Psi2S_Px, fit_Psi2S_Py, fit_Psi2S_Pz, fit_Psi2S_M, fit_Psi2S_ChiSq;
332  std::vector<float> fit_Upsi1S_Px, fit_Upsi1S_Py, fit_Upsi1S_Pz, fit_Upsi1S_M, fit_Upsi1S_ChiSq;
333  std::vector<float> fit_Upsi2S_Px, fit_Upsi2S_Py, fit_Upsi2S_Pz, fit_Upsi2S_M, fit_Upsi2S_ChiSq;
334  std::vector<float> fit_Upsi3S_Px, fit_Upsi3S_Py, fit_Upsi3S_Pz, fit_Upsi3S_M, fit_Upsi3S_ChiSq;
335 
336  // Loop over di-muon vertices associated with a candidate
337  for(std::vector<const xAOD::Vertex*>::const_iterator vtxItr = candidateOniaVertices.begin(); vtxItr != candidateOniaVertices.end(); ++vtxItr ) {
338 
339  //------------------------------------
340  // Add an element link to each dimuon which formed a
341  // candidate, leading to the decision to save this conversion
342  //------------------------------------
343  ElementLink<xAOD::VertexContainer> myLink;
344  myLink.setElement(*vtxItr);
345  myLink.setStorableObject(*diMuonContainer);
346 
347  if(!myLink.isValid()) {
348  ATH_MSG_WARNING("Invalid DiMuon ElementLink!");
349  }
350 
351  diMuonLinks.push_back(myLink);
352  //------------------------------------
353 
354  // Check which mass window this di-muon passed
355  static const SG::ConstAccessor<Char_t> passed_PsiAcc("passed_Psi");
356  static const SG::ConstAccessor<Char_t> passed_UpsiAcc("passed_Upsi");
357  bool passed_Psi = passed_PsiAcc(**vtxItr);
358  bool passed_Upsi = passed_UpsiAcc(**vtxItr);
359 
360  //------------------------------------
361  // Cascade fit with J/psi mass hypothesis
362  //------------------------------------
363  float fitChiSq_Psi1S = 99999;
364  TLorentzVector fitResult_Psi1S;
365 
366  // Only bother with the fit if di-muon mass is within the relveant range,
367  // but still fill an dummy 4-vector to preserve one to one correspondance with "DiMuonLinks"
368  if(passed_Psi) {
369  ATH_CHECK( doCascadeFit(*vtxItr,constConvVertex,3096.916,fitResult_Psi1S,fitChiSq_Psi1S) );
370  }
371 
372  fit_Psi1S_Px.push_back(fitResult_Psi1S.Px());
373  fit_Psi1S_Py.push_back(fitResult_Psi1S.Py());
374  fit_Psi1S_Pz.push_back(fitResult_Psi1S.Pz());
375  fit_Psi1S_M.push_back(fitResult_Psi1S.M());
376  fit_Psi1S_ChiSq.push_back(fitChiSq_Psi1S);
377 
378  //------------------------------------
379  // Cascade fit with psi(2S) mass hypothesis
380  //------------------------------------
381  float fitChiSq_Psi2S = 99999;
382  TLorentzVector fitResult_Psi2S;
383 
384  // Only bother with the fit if di-muon mass is within the relveant range,
385  // but still fill an dummy 4-vector to preserve one to one correspondance with "DiMuonLinks"
386  if(passed_Psi) {
387  ATH_CHECK( doCascadeFit(*vtxItr,constConvVertex,3686.097,fitResult_Psi2S,fitChiSq_Psi2S) );
388  }
389 
390  fit_Psi2S_Px.push_back(fitResult_Psi2S.Px());
391  fit_Psi2S_Py.push_back(fitResult_Psi2S.Py());
392  fit_Psi2S_Pz.push_back(fitResult_Psi2S.Pz());
393  fit_Psi2S_M.push_back(fitResult_Psi2S.M());
394  fit_Psi2S_ChiSq.push_back(fitChiSq_Psi2S);
395 
396  //------------------------------------
397  // Cascade fit with Upsi(1S) mass hypothesis
398  //------------------------------------
399  float fitChiSq_Upsi1S = 99999;
400  TLorentzVector fitResult_Upsi1S;
401 
402  // Only bother with the fit if di-muon mass is within the relveant range,
403  // but still fill an dummy 4-vector to preserve one to one correspondance with "DiMuonLinks"
404  if(passed_Upsi) {
405  ATH_CHECK( doCascadeFit(*vtxItr,constConvVertex,9460.30,fitResult_Upsi1S,fitChiSq_Upsi1S) );
406  }
407 
408  fit_Upsi1S_Px.push_back(fitResult_Upsi1S.Px());
409  fit_Upsi1S_Py.push_back(fitResult_Upsi1S.Py());
410  fit_Upsi1S_Pz.push_back(fitResult_Upsi1S.Pz());
411  fit_Upsi1S_M.push_back(fitResult_Upsi1S.M());
412  fit_Upsi1S_ChiSq.push_back(fitChiSq_Upsi1S);
413 
414  //------------------------------------
415  // Cascade fit with Upsi(2S) mass hypothesis
416  //------------------------------------
417  float fitChiSq_Upsi2S = 99999;
418  TLorentzVector fitResult_Upsi2S;
419 
420  // Only bother with the fit if di-muon mass is within the relveant range,
421  // but still fill an dummy 4-vector to preserve one to one correspondance with "DiMuonLinks"
422  if(passed_Upsi) {
423  ATH_CHECK( doCascadeFit(*vtxItr,constConvVertex,10023.26,fitResult_Upsi2S,fitChiSq_Upsi2S) );
424  }
425 
426  fit_Upsi2S_Px.push_back(fitResult_Upsi2S.Px());
427  fit_Upsi2S_Py.push_back(fitResult_Upsi2S.Py());
428  fit_Upsi2S_Pz.push_back(fitResult_Upsi2S.Pz());
429  fit_Upsi2S_M.push_back(fitResult_Upsi2S.M());
430  fit_Upsi2S_ChiSq.push_back(fitChiSq_Upsi2S);
431 
432  //------------------------------------
433  // Cascade fit with Upsi(3S) mass hypothesis
434  //------------------------------------
435  float fitChiSq_Upsi3S = 99999;
436  TLorentzVector fitResult_Upsi3S;
437 
438  // Only bother with the fit if di-muon mass is within the relveant range,
439  // but still fill an dummy 4-vector to preserve one to one correspondance with "DiMuonLinks"
440  if(passed_Upsi) {
441  ATH_CHECK( doCascadeFit(*vtxItr,constConvVertex,10355.2,fitResult_Upsi3S,fitChiSq_Upsi3S) );
442  }
443 
444  fit_Upsi3S_Px.push_back(fitResult_Upsi3S.Px());
445  fit_Upsi3S_Py.push_back(fitResult_Upsi3S.Py());
446  fit_Upsi3S_Pz.push_back(fitResult_Upsi3S.Pz());
447  fit_Upsi3S_M.push_back(fitResult_Upsi3S.M());
448  fit_Upsi3S_ChiSq.push_back(fitChiSq_Upsi3S);
449 
450  }
451 
452  //------------------------------------
453  // Decorate selected conversions
454  //------------------------------------
455  ATH_MSG_DEBUG("Decorating conversion vertices");
456 
457  static const SG::Accessor< std::vector< ElementLink<xAOD::VertexContainer> > > DiMuonLinksAcc("DiMuonLinks");
458  DiMuonLinksAcc(*convVertexCandidate) = diMuonLinks;
459 
460  static const SG::Accessor< std::vector<float> > CascadeFit_Psi1S_PxAcc("CascadeFit_Psi1S_Px");
461  static const SG::Accessor< std::vector<float> > CascadeFit_Psi1S_PyAcc("CascadeFit_Psi1S_Py");
462  static const SG::Accessor< std::vector<float> > CascadeFit_Psi1S_PzAcc("CascadeFit_Psi1S_Pz");
463  static const SG::Accessor< std::vector<float> > CascadeFit_Psi1S_MAcc("CascadeFit_Psi1S_M");
464  static const SG::Accessor< std::vector<float> > CascadeFit_Psi1S_ChiSqAcc("CascadeFit_Psi1S_ChiSq");
465  CascadeFit_Psi1S_PxAcc(*convVertexCandidate) = fit_Psi1S_Px;
466  CascadeFit_Psi1S_PyAcc(*convVertexCandidate) = fit_Psi1S_Py;
467  CascadeFit_Psi1S_PzAcc(*convVertexCandidate) = fit_Psi1S_Pz;
468  CascadeFit_Psi1S_MAcc(*convVertexCandidate) = fit_Psi1S_M;
469  CascadeFit_Psi1S_ChiSqAcc(*convVertexCandidate) = fit_Psi1S_ChiSq;
470 
471  static const SG::Accessor< std::vector<float> > CascadeFit_Psi2S_PxAcc("CascadeFit_Psi2S_Px");
472  static const SG::Accessor< std::vector<float> > CascadeFit_Psi2S_PyAcc("CascadeFit_Psi2S_Py");
473  static const SG::Accessor< std::vector<float> > CascadeFit_Psi2S_PzAcc("CascadeFit_Psi2S_Pz");
474  static const SG::Accessor< std::vector<float> > CascadeFit_Psi2S_MAcc("CascadeFit_Psi2S_M");
475  static const SG::Accessor< std::vector<float> > CascadeFit_Psi2S_ChiSqAcc("CascadeFit_Psi2S_ChiSq");
476  CascadeFit_Psi2S_PxAcc(*convVertexCandidate) = fit_Psi2S_Px;
477  CascadeFit_Psi2S_PyAcc(*convVertexCandidate) = fit_Psi2S_Py;
478  CascadeFit_Psi2S_PzAcc(*convVertexCandidate) = fit_Psi2S_Pz;
479  CascadeFit_Psi2S_MAcc(*convVertexCandidate) = fit_Psi2S_M;
480  CascadeFit_Psi2S_ChiSqAcc(*convVertexCandidate) = fit_Psi2S_ChiSq;
481 
482  static const SG::Accessor< std::vector<float> > CascadeFit_Upsi1S_PxAcc("CascadeFit_Upsi1S_Px");
483  static const SG::Accessor< std::vector<float> > CascadeFit_Upsi1S_PyAcc("CascadeFit_Upsi1S_Py");
484  static const SG::Accessor< std::vector<float> > CascadeFit_Upsi1S_PzAcc("CascadeFit_Upsi1S_Pz");
485  static const SG::Accessor< std::vector<float> > CascadeFit_Upsi1S_MAcc("CascadeFit_Upsi1S_M");
486  static const SG::Accessor< std::vector<float> > CascadeFit_Upsi1S_ChiSqAcc("CascadeFit_Upsi1S_ChiSq");
487  CascadeFit_Upsi1S_PxAcc(*convVertexCandidate) = fit_Upsi1S_Px;
488  CascadeFit_Upsi1S_PyAcc(*convVertexCandidate) = fit_Upsi1S_Py;
489  CascadeFit_Upsi1S_PzAcc(*convVertexCandidate) = fit_Upsi1S_Pz;
490  CascadeFit_Upsi1S_MAcc(*convVertexCandidate) = fit_Upsi1S_M;
491  CascadeFit_Upsi1S_ChiSqAcc(*convVertexCandidate) = fit_Upsi1S_ChiSq;
492 
493  static const SG::Accessor< std::vector<float> > CascadeFit_Upsi2S_PxAcc("CascadeFit_Upsi2S_Px");
494  static const SG::Accessor< std::vector<float> > CascadeFit_Upsi2S_PyAcc("CascadeFit_Upsi2S_Py");
495  static const SG::Accessor< std::vector<float> > CascadeFit_Upsi2S_PzAcc("CascadeFit_Upsi2S_Pz");
496  static const SG::Accessor< std::vector<float> > CascadeFit_Upsi2S_MAcc("CascadeFit_Upsi2S_M");
497  static const SG::Accessor< std::vector<float> > CascadeFit_Upsi2S_ChiSqAcc("CascadeFit_Upsi2S_ChiSq");
498  CascadeFit_Upsi2S_PxAcc(*convVertexCandidate) = fit_Upsi2S_Px;
499  CascadeFit_Upsi2S_PyAcc(*convVertexCandidate) = fit_Upsi2S_Py;
500  CascadeFit_Upsi2S_PzAcc(*convVertexCandidate) = fit_Upsi2S_Pz;
501  CascadeFit_Upsi2S_MAcc(*convVertexCandidate) = fit_Upsi2S_M;
502  CascadeFit_Upsi2S_ChiSqAcc(*convVertexCandidate) = fit_Upsi2S_ChiSq;
503 
504  static const SG::Accessor< std::vector<float> > CascadeFit_Upsi3S_PxAcc("CascadeFit_Upsi3S_Px");
505  static const SG::Accessor< std::vector<float> > CascadeFit_Upsi3S_PyAcc("CascadeFit_Upsi3S_Py");
506  static const SG::Accessor< std::vector<float> > CascadeFit_Upsi3S_PzAcc("CascadeFit_Upsi3S_Pz");
507  static const SG::Accessor< std::vector<float> > CascadeFit_Upsi3S_MAcc("CascadeFit_Upsi3S_M");
508  static const SG::Accessor< std::vector<float> > CascadeFit_Upsi3S_ChiSqAcc("CascadeFit_Upsi3S_ChiSq");
509  CascadeFit_Upsi3S_PxAcc(*convVertexCandidate) = fit_Upsi3S_Px;
510  CascadeFit_Upsi3S_PyAcc(*convVertexCandidate) = fit_Upsi3S_Py;
511  CascadeFit_Upsi3S_PzAcc(*convVertexCandidate) = fit_Upsi3S_Pz;
512  CascadeFit_Upsi3S_MAcc(*convVertexCandidate) = fit_Upsi3S_M;
513  CascadeFit_Upsi3S_ChiSqAcc(*convVertexCandidate) = fit_Upsi3S_ChiSq;
514 
515  static const SG::Accessor<float> pxAcc("px");
516  static const SG::Accessor<float> pyAcc("py");
517  static const SG::Accessor<float> pzAcc("pz");
518  pxAcc(*convVertexCandidate) = momentum.x();
519  pyAcc(*convVertexCandidate) = momentum.y();
520  pzAcc(*convVertexCandidate) = momentum.z();
521 
522  static const SG::Accessor<float> deltaCotThetaTrkAcc("deltaCotThetaTrk");
523  static const SG::Accessor<float> minimumDistanceTrkAcc("minimumDistanceTrk");
524  deltaCotThetaTrkAcc(*convVertexCandidate) = deltaCotTheta;
525  minimumDistanceTrkAcc(*convVertexCandidate) = distance;
526 
527  static const SG::Accessor<float> deltaPhiTracksAcc("deltaPhiTracks");
528  static const SG::Accessor<float> DR1R2Acc("DR1R2");
529  deltaPhiTracksAcc(*convVertexCandidate) = vertexOutput["deltaPhiTracks"];
530  DR1R2Acc(*convVertexCandidate) = vertexOutput["DR1R2"];
531 
532  static const SG::Accessor<Char_t> passedAcc("passed");
533  passedAcc(*convVertexCandidate) = true; // Used in event skimming
534 
535  conversionContainer->push_back(convVertexCandidate.release());
536 
537  } else {
538  ATH_MSG_DEBUG("Vertex Fit Failed");
539  }
540 
541  } // Neg Track Loop
542 
543  } // Pos Track Loop
544 
545  } // callConvFinder
546 
547  // Write the results to StoreGate
548  CHECK(evtStore()->record(conversionContainer.release(), m_conversionContainerName));
549  CHECK(evtStore()->record(conversionAuxContainer.release(), m_conversionContainerName+"Aux."));
550 
551  ATH_MSG_DEBUG("-------------------------");
552  ATH_MSG_DEBUG("Number of track pairs: " << nTrackPairs_Init);
553  ATH_MSG_DEBUG("Number of track pairs selected: " << nTrackPairs_Selected);
554  ATH_MSG_DEBUG("Number of successful vertex fits: " << nConv_VertexFit);
555  ATH_MSG_DEBUG("Number of selected vertices: " << nConv_Selected);
556  ATH_MSG_DEBUG("Number of selected vertices (after DeltaM req.): " << nConv_Selected_DeltaM);
557 
558  return StatusCode::SUCCESS;
559  }
560 
561  StatusCode BPhysConversionFinder::doCascadeFit(const xAOD::Vertex * diMuonVertex, const xAOD::Vertex * convVertex, const double diMuonMassConstraint, TLorentzVector & fitMom, float & chiSq) const {
562 
563  std::vector<const xAOD::TrackParticle*> diMuonTracks;
564  diMuonTracks.push_back(diMuonVertex->trackParticle(0));
565  diMuonTracks.push_back(diMuonVertex->trackParticle(1));
566 
567  std::vector<double> diMuonTrackMasses;
568  diMuonTrackMasses.push_back(105.658);
569  diMuonTrackMasses.push_back(105.658);
570 
571  std::vector<const xAOD::TrackParticle*> convTracks;
572  convTracks.push_back(convVertex->trackParticle(0));
573  convTracks.push_back(convVertex->trackParticle(1));
574 
575  std::vector<double> convTrackMasses;
576  convTrackMasses.push_back(0.511);
577  convTrackMasses.push_back(0.511);
578 
579  // Reset
580  std::unique_ptr<Trk::IVKalState> state = m_cascadeFitter->makeState();
581 
582  // Set Robustness
583  m_cascadeFitter->setRobustness(0, *state);
584 
585  // Build up the topology
586 
587  // Vertex list
588  std::vector<Trk::VertexID> vrtList;
589  // V0 vertex
590  Trk::VertexID vID;
591  vID = m_cascadeFitter->startVertex(convTracks,convTrackMasses,*state,0.0); // Constrain converision mass to zero
592 
593  vrtList.push_back(vID);
594 
595  // chi_c/b vertex
596  Trk::VertexID vID2 = m_cascadeFitter->nextVertex(diMuonTracks,diMuonTrackMasses,vrtList,*state);
597 
598  std::vector<Trk::VertexID> cnstV;
599  cnstV.clear();
600  if ( !m_cascadeFitter->addMassConstraint(vID2,diMuonTracks,cnstV,*state,diMuonMassConstraint).isSuccess() ) {
601  ATH_MSG_WARNING("addMassConstraint failed");
602  }
603 
604  // Do the fit
605  std::unique_ptr<Trk::VxCascadeInfo> result(m_cascadeFitter->fitCascade(*state));
606 
607  const std::vector< std::vector<TLorentzVector> > &moms = result->getParticleMoms();
608 
609  {
610 
611  if(moms.size() > 2) ATH_MSG_WARNING("DoCascadeFit - More than two output momentum!?");
612 
613  TLorentzVector conv_Fit = moms.at(0).at(0) + moms.at(0).at(1);
614  TLorentzVector diMuon_Fit = moms.at(1).at(0) + moms.at(1).at(1);
615 
616  // Momentum of DiMuon + photon system
617  fitMom = diMuon_Fit + conv_Fit;
618 
619  chiSq = result->fitChi2()/result->nDoF();
620 
621  // Done with the fit result
622  result.reset();
623 
624  }
625 
626  return StatusCode::SUCCESS;
627 
628  }
629 
630 
631 }
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
DerivationFramework::BPhysConversionFinder::m_v0Tools
ToolHandle< Trk::V0Tools > m_v0Tools
Definition: BPhysConversionFinder.h:51
V0Tools.h
DataModel_detail::const_iterator
Const iterator class for DataVector/DataList.
Definition: DVLIterator.h:82
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
Trk::VertexID
int VertexID
Definition: IVertexCascadeFitter.h:23
xAOD::uint8_t
uint8_t
Definition: Muon_v1.cxx:575
DerivationFramework::BPhysConversionFinder::m_conversionContainerName
std::string m_conversionContainerName
Definition: BPhysConversionFinder.h:59
SG::Accessor
Helper class to provide type-safe access to aux data.
Definition: Control/AthContainers/AthContainers/Accessor.h:66
DerivationFramework::BPhysConversionFinder::m_maxDistBetweenTracks
float m_maxDistBetweenTracks
Definition: BPhysConversionFinder.h:61
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
Trk::ParametersT
Dummy class used to allow special convertors to be called for surfaces owned by a detector element.
Definition: EMErrorDetail.h:25
DerivationFramework::BPhysConversionFinder::m_diMuonCollectionToCheck
std::string m_diMuonCollectionToCheck
Definition: BPhysConversionFinder.h:48
xAOD::numberOfPixelHits
@ numberOfPixelHits
these are the pixel hits, including the b-layer [unit8_t].
Definition: TrackingPrimitives.h:259
SG::ConstAccessor
Helper class to provide constant type-safe access to aux data.
Definition: ConstAccessor.h:54
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
DerivationFramework::BPhysConversionFinder::doCascadeFit
StatusCode doCascadeFit(const xAOD::Vertex *diMuonVertex, const xAOD::Vertex *convVertex, const double diMuonMassConstraint, TLorentzVector &fitMom, float &chiSq) const
Definition: BPhysConversionFinder.cxx:561
DerivationFramework::BPhysConversionFinder::m_distanceTool
ToolHandle< Trk::ITrkDistanceFinder > m_distanceTool
Definition: BPhysConversionFinder.h:54
TrkVKalVrtFitter.h
DerivationFramework::BPhysConversionFinder::m_cascadeFitter
ToolHandle< Trk::TrkVKalVrtFitter > m_cascadeFitter
Definition: BPhysConversionFinder.h:56
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::Vertex_v1::addTrackAtVertex
void addTrackAtVertex(const ElementLink< TrackParticleContainer > &tr, float weight=1.0)
Add a new track to the vertex.
Definition: Vertex_v1.cxx:314
xAOD::TrackParticle_v1::perigeeParameters
const Trk::Perigee & perigeeParameters() const
Returns the Trk::MeasuredPerigee track parameters.
Definition: TrackParticle_v1.cxx:485
DerivationFramework::BPhysConversionFinder::m_vertexFitter
ToolHandle< Trk::IVertexFitter > m_vertexFitter
Definition: BPhysConversionFinder.h:52
DerivationFramework::BPhysConversionFinder::m_passFlagsToCheck
std::vector< std::string > m_passFlagsToCheck
Definition: BPhysConversionFinder.h:49
ParticleGun_EoverP_Config.momentum
momentum
Definition: ParticleGun_EoverP_Config.py:63
DerivationFramework::BPhysConversionFinder::initialize
StatusCode initialize() override
Definition: BPhysConversionFinder.cxx:61
beamspotman.n
n
Definition: beamspotman.py:731
Trk::theta
@ theta
Definition: ParamDefs.h:72
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
master.flag
bool flag
Definition: master.py:29
AthenaPoolTestRead.acc
acc
Definition: AthenaPoolTestRead.py:16
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
CHECK
#define CHECK(...)
Evaluate an expression and check for errors.
Definition: Control/AthenaKernel/AthenaKernel/errorcheck.h:422
DerivationFramework::BPhysConversionFinder::BPhysConversionFinder
BPhysConversionFinder(const std::string &t, const std::string &n, const IInterface *p)
Definition: BPhysConversionFinder.cxx:22
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
xAOD::Vertex_v1::trackParticle
const TrackParticle * trackParticle(size_t i) const
Get the pointer to a given track that was used in vertex reco.
Definition: Vertex_v1.cxx:249
xAOD::Vertex_v1::clearTracks
void clearTracks()
Remove all tracks from the vertex.
Definition: Vertex_v1.cxx:331
DerivationFramework
THE reconstruction tool.
Definition: ParticleSortingAlg.h:24
SG::AuxElement::index
size_t index() const
Return the index of this element within its container.
DataVector
Derived DataVector<T>.
Definition: DataVector.h:581
DerivationFramework::BPhysConversionFinder::m_requireDeltaM
bool m_requireDeltaM
Definition: BPhysConversionFinder.h:64
DerivationFramework::BPhysConversionFinder::m_vertexEstimator
ToolHandle< InDet::VertexPointEstimator > m_vertexEstimator
Definition: BPhysConversionFinder.h:53
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
IVertexFitter.h
DerivationFramework::BPhysConversionFinder::m_postSelector
ToolHandle< InDet::ConversionPostSelector > m_postSelector
Definition: BPhysConversionFinder.h:55
VxCascadeInfo.h
DataVector::end
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
Trk::vertex
@ vertex
Definition: MeasurementType.h:21
VertexContainer.h
xAOD::photon
@ photon
Definition: TrackingPrimitives.h:199
xAOD::Vertex_v1
Class describing a Vertex.
Definition: Vertex_v1.h:42
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
GeoPrimitivesHelpers.h
DerivationFramework::BPhysConversionFinder::addBranches
virtual StatusCode addBranches() const override
Pass the thinning service
Definition: BPhysConversionFinder.cxx:87
DerivationFramework::BPhysConversionFinder::m_inputTrackParticleContainerName
std::string m_inputTrackParticleContainerName
Definition: BPhysConversionFinder.h:58
xAOD::numberOfSCTHits
@ numberOfSCTHits
number of hits in SCT [unit8_t].
Definition: TrackingPrimitives.h:268
DerivationFramework::BPhysConversionFinder::m_maxDeltaCotTheta
float m_maxDeltaCotTheta
Definition: BPhysConversionFinder.h:62
DerivationFramework::BPhysConversionFinder::m_maxDeltaM
float m_maxDeltaM
Definition: BPhysConversionFinder.h:65
xAOD::track
@ track
Definition: TrackingPrimitives.h:512
xAOD::TrackParticle_v1
Class describing a TrackParticle.
Definition: TrackParticle_v1.h:43
ConstAccessor.h
Helper class to provide constant type-safe access to aux data.
DerivationFramework::BPhysConversionFinder::finalize
StatusCode finalize() override
Definition: BPhysConversionFinder.cxx:79
AthAlgTool
Definition: AthAlgTool.h:26
BPhysConversionFinder.h
Amg::distance
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
Definition: GeoPrimitivesHelpers.h:54
DataVector::size
size_type size() const noexcept
Returns the number of elements in the collection.
DataVector::begin
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
VertexAuxContainer.h
Generated on Tue Jul 2 2024 21:11:06 for ATLAS Offline Software by   doxygen 1.8.18