ATLAS Offline Software
Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
7 // Header include
8 #include "NewVrtSecInclusiveTool/NewVrtSecInclusiveTool.h"
9 #include "TrkNeutralParameters/NeutralParameters.h"
10 #include "TrkTrackSummary/TrackSummary.h"
11 #include "TrkVKalVrtFitter/TrkVKalVrtFitter.h"
12 #include "AthContainers/ConstAccessor.h"
13 #include "CxxUtils/sincos.h"
14 //-------------------------------------------------
15 #include "TrkGeometry/TrackingGeometry.h"
16 #include "TrkGeometry/TrackingVolume.h"
17 #include "TrkGeometry/Layer.h"
18 #include "TruthUtils/MagicNumbers.h"
19 #include "GeoPrimitives/GeoPrimitivesHelpers.h"
20  // Other stuff
21 #include <cmath>
22 
23 
24 namespace Rec{
25 
26 
27  void NewVrtSecInclusiveTool::printWrkSet(const std::vector<WrkVrt> *WrkVrtSet, const std::string &name) const {
28  int nGoodV=0;
29  if(msgLvl(MSG::INFO)){
30  for(int iv=0; iv<(int)WrkVrtSet->size(); iv++) {
31  msg(MSG::INFO)<<name
32  <<"= "<<(*WrkVrtSet)[iv].vertex[0]
33  <<", "<<(*WrkVrtSet)[iv].vertex[1]
34  <<", "<<(*WrkVrtSet)[iv].vertex[2]
35  <<" NTrk="<<(*WrkVrtSet)[iv].selTrk.size()
36  <<" is good="<<std::boolalpha<<(*WrkVrtSet)[iv].Good<<std::noboolalpha
37  <<" Chi2="<<(*WrkVrtSet)[iv].chi2
38  <<" Mass="<<(*WrkVrtSet)[iv].vertexMom.M()
39  <<" detached="<<(*WrkVrtSet)[iv].detachedTrack
40  <<" proj.dist="<<(*WrkVrtSet)[iv].projectedVrt
41  <<" trk=";
42  for(int kk=0; kk<(int)(*WrkVrtSet)[iv].selTrk.size(); kk++) {
43  msg(MSG::INFO)<<", "<<(*WrkVrtSet)[iv].selTrk[kk];}
44  for(int kk=0; kk<(int)(*WrkVrtSet)[iv].selTrk.size(); kk++) {
45  msg(MSG::INFO)<<", "<<momAtVrt((*WrkVrtSet)[iv].trkAtVrt[kk]).Pt();}
46  msg(MSG::INFO)<<endmsg;
47  if((*WrkVrtSet)[iv].Good)nGoodV++;
48  }
49  msg(MSG::INFO)<<name<<" N="<<nGoodV<<endmsg;
50  }
51  }
52 
53  /* Technicalities */
54  double NewVrtSecInclusiveTool::projSV_PV(const Amg::Vector3D & SV, const xAOD::Vertex & PV, const TLorentzVector & Direction)
55  {
56  TVector3 SV_PV( SV.x()-PV.x(), SV.y()-PV.y(), SV.z()-PV.z() );
57  return Direction.Vect().Unit()*SV_PV.Unit();
58  }
59 
60  double NewVrtSecInclusiveTool::vrtVrtDist(const xAOD::Vertex & primVrt, const Amg::Vector3D & secVrt,
61  const std::vector<double>& secVrtErr, double& signif)
62 
63  {
64  double distx = primVrt.x()- secVrt.x();
65  double disty = primVrt.y()- secVrt.y();
66  double distz = primVrt.z()- secVrt.z();
67 
68 
69  AmgSymMatrix(3) primCovMtx=primVrt.covariancePosition(); //Create
70  primCovMtx(0,0) += secVrtErr[0];
71  primCovMtx(0,1) += secVrtErr[1];
72  primCovMtx(1,0) += secVrtErr[1];
73  primCovMtx(1,1) += secVrtErr[2];
74  primCovMtx(0,2) += secVrtErr[3];
75  primCovMtx(2,0) += secVrtErr[3];
76  primCovMtx(1,2) += secVrtErr[4];
77  primCovMtx(2,1) += secVrtErr[4];
78  primCovMtx(2,2) += secVrtErr[5];
79 
80  AmgSymMatrix(3) wgtMtx = primCovMtx.inverse();
81 
82  signif = distx*wgtMtx(0,0)*distx
83  +disty*wgtMtx(1,1)*disty
84  +distz*wgtMtx(2,2)*distz
85  +2.*distx*wgtMtx(0,1)*disty
86  +2.*distx*wgtMtx(0,2)*distz
87  +2.*disty*wgtMtx(1,2)*distz;
88  signif=std::sqrt(std::abs(signif));
89  if( signif!=signif ) signif = 0.;
90  return std::sqrt(distx*distx+disty*disty+distz*distz);
91  }
92 
93  double NewVrtSecInclusiveTool::vrtVrtDist2D(const xAOD::Vertex & primVrt, const Amg::Vector3D & secVrt,
94  const std::vector<double>& secVrtErr, double& signif)
95 
96  {
97  double distx = primVrt.x()- secVrt.x();
98  double disty = primVrt.y()- secVrt.y();
99 
100 
101  AmgSymMatrix(3) primCovMtx=primVrt.covariancePosition(); //Create
102  AmgSymMatrix(2) covMtx;
103  covMtx(0,0) = primCovMtx(0,0) + secVrtErr[0];
104  covMtx(0,1) = primCovMtx(0,1) + secVrtErr[1];
105  covMtx(1,0) = primCovMtx(1,0) + secVrtErr[1];
106  covMtx(1,1) = primCovMtx(1,1) + secVrtErr[2];
107 
108  AmgSymMatrix(2) wgtMtx = covMtx.inverse();
109 
110  signif = distx*wgtMtx(0,0)*distx
111  +disty*wgtMtx(1,1)*disty
112  +2.*distx*wgtMtx(0,1)*disty;
113  signif=std::sqrt(std::abs(signif));
114  if( signif!=signif ) signif = 0.;
115  return std::sqrt(distx*distx+disty*disty);
116  }
117 
118 
119  double NewVrtSecInclusiveTool::vrtVrtDist(const Amg::Vector3D & vrt1, const std::vector<double> & vrtErr1,
120  const Amg::Vector3D & vrt2, const std::vector<double> & vrtErr2)
121 
122  {
123  double distx = vrt1.x()- vrt2.x();
124  double disty = vrt1.y()- vrt2.y();
125  double distz = vrt1.z()- vrt2.z();
126 
127  AmgSymMatrix(3) primCovMtx; //Create
128  primCovMtx(0,0) = vrtErr1[0]+vrtErr2[0];
129  primCovMtx(0,1) = primCovMtx(1,0) = vrtErr1[1]+vrtErr2[1];
130  primCovMtx(1,1) = vrtErr1[2]+vrtErr2[2];
131  primCovMtx(0,2) = primCovMtx(2,0) = vrtErr1[3]+vrtErr2[3];
132  primCovMtx(1,2) = primCovMtx(2,1) = vrtErr1[4]+vrtErr2[4];
133  primCovMtx(2,2) = vrtErr1[5]+vrtErr2[5];
134 
135  AmgSymMatrix(3) wgtMtx = primCovMtx.inverse();
136 
137  double signif =
138  distx*wgtMtx(0,0)*distx
139  +disty*wgtMtx(1,1)*disty
140  +distz*wgtMtx(2,2)*distz
141  +2.*distx*wgtMtx(0,1)*disty
142  +2.*distx*wgtMtx(0,2)*distz
143  +2.*disty*wgtMtx(1,2)*distz;
144  signif=std::sqrt(std::abs(signif));
145  if(signif != signif) signif = 0.;
146  return signif;
147  }
148 //
149  double NewVrtSecInclusiveTool::PntPntDist(const Amg::Vector3D & Vrt1, const Amg::Vector3D & Vrt2)
150  {
151  double dx = Vrt1.x()- Vrt2.x();
152  double dy = Vrt1.y()- Vrt2.y();
153  double dz = Vrt1.z()- Vrt2.z();
154  return std::sqrt(dx*dx+dy*dy*dz*dz);
155  }
156 
157 //----------------------------
158 // Vertex error along radius
159 //----------------------------
160  double NewVrtSecInclusiveTool::vrtRadiusError(const Amg::Vector3D & SecVrt, const std::vector<double> & VrtErr)
161  {
162  double DirX=SecVrt.x(), DirY=SecVrt.y();
163  double Covar = DirX*VrtErr[0]*DirX
164  +2.*DirX*VrtErr[1]*DirY
165  +DirY*VrtErr[2]*DirY;
166  Covar /= DirX*DirX + DirY*DirY;
167  Covar=std::sqrt(std::abs(Covar));
168  if(Covar != Covar) Covar = 0.;
169  return Covar;
170  }
171 
172 
173  /* Invariant mass calculation for V0 decays*/
174  /* Gives correct mass assignment in case of nonequal masses*/
175 
176 
177  double NewVrtSecInclusiveTool::massV0(const std::vector< std::vector<double> >& TrkAtVrt,
178  double massP, double massPi )
179 
180  {
181  double ap1i=std::abs(TrkAtVrt[0][2]); double ap2i=std::abs(TrkAtVrt[1][2]);
182  CxxUtils::sincos phi1(TrkAtVrt[0][0]);
183  CxxUtils::sincos theta1(TrkAtVrt[0][1]);
184  CxxUtils::sincos phi2(TrkAtVrt[1][0]);
185  CxxUtils::sincos theta2(TrkAtVrt[1][1]);
186  double px = phi1.cs * theta1.sn * ap1i
187  + phi2.cs * theta2.sn * ap2i;
188  double py = phi1.sn * theta1.sn * ap1i
189  + phi2.sn * theta2.sn * ap2i;
190  double pz = theta1.cs * ap1i
191  + theta2.cs * ap2i;
192  double ee= (ap1i > ap2i) ?
193  (std::sqrt(ap1i*ap1i+massP*massP)+std::sqrt(ap2i*ap2i+massPi*massPi)):
194  (std::sqrt(ap2i*ap2i+massP*massP)+std::sqrt(ap1i*ap1i+massPi*massPi));
195  double test=(ee-pz)*(ee+pz)-px*px-py*py;
196  return test>0 ? std::sqrt(test) : 0.;
197  }
198 
199 
200 
201  ROOT::Math::PxPyPzEVector NewVrtSecInclusiveTool::momAtVrt(const std::vector< double >& inpTrk) const
202  {
203  double api=1./std::abs(inpTrk[2]);
204  CxxUtils::sincos phi(inpTrk[0]);
205  CxxUtils::sincos theta(inpTrk[1]);
206  double px = phi.cs * theta.sn * api;
207  double py = phi.sn * theta.sn * api;
208  double pz = theta.cs * api;
209  double ee = std::sqrt( px*px + py*py + pz*pz + m_massPi*m_massPi);
210  return {px,py,pz,ee};
211  }
212 
213 /*************************************************************************************************************/
214  int NewVrtSecInclusiveTool::getIBLHit(const xAOD::TrackParticle* Part)
215  {
216  uint8_t IBLhit,IBLexp;
217  if(!Part->summaryValue( IBLexp, xAOD::expectInnermostPixelLayerHit) ) IBLexp = 0;
218  if( IBLexp==0 ) return -1;
219  if(!Part->summaryValue( IBLhit, xAOD::numberOfInnermostPixelLayerHits) ) IBLhit = 0;
220  if(IBLhit) return 1;
221  else return 0;
222  }
223  int NewVrtSecInclusiveTool::getBLHit(const xAOD::TrackParticle* Part)
224  {
225  uint8_t BLhit,BLexp;
226  if(!Part->summaryValue( BLexp, xAOD::expectNextToInnermostPixelLayerHit) ) BLexp = 0;
227  if( BLexp==0 ) return -1;
228  if(!Part->summaryValue( BLhit, xAOD::numberOfNextToInnermostPixelLayerHits) ) BLhit = 0;
229  if(BLhit) return 1;
230  else return 0;
231  }
232 
233  void NewVrtSecInclusiveTool::getPixelDiscs(const xAOD::TrackParticle* Part, int &d0Hit, int &d1Hit, int &d2Hit)
234  {
235  uint32_t HitPattern=Part->hitPattern();
236  d0Hit=0; if( HitPattern&((1<<Trk::pixelEndCap0)) ) d0Hit=1;
237  d1Hit=0; if( HitPattern&((1<<Trk::pixelEndCap1)) ) d1Hit=1;
238  d2Hit=0; if( HitPattern&((1<<Trk::pixelEndCap2)) ) d2Hit=1;
239  }
240 /*************************************************************************************************************/
241 
242  int NewVrtSecInclusiveTool::getIdHF(const xAOD::TrackParticle* TP ) const {
243  static const SG::ConstAccessor< ElementLink< xAOD::TruthParticleContainer> >
244  truthParticleLinkAcc ( "truthParticleLink");
245  if( truthParticleLinkAcc.isAvailable(*TP) ) {
246  const ElementLink<xAOD::TruthParticleContainer>& tplink = truthParticleLinkAcc(*TP);
247  if( !tplink.isValid() ) return 0;
248  static const SG::ConstAccessor< float >truthMatchProbabilityAcc ( "truthMatchProbability" );
249  if( truthMatchProbabilityAcc( *TP ) < 0.75 ) return 0;
250  if( HepMC::is_simulation_particle((*tplink))) return 0;
251  //-------------
252  const xAOD::TruthParticle * parent1=getPreviousParent(*tplink);
253  if(!parent1 || !parent1->isHadron()) return 0;
254  if(parent1->isBottomHadron()) return isExcitedHadron(parent1) ? 0 : 1; //Ground state B-hadron
255  if(parent1->isCharmHadron() && !isExcitedHadron(parent1))return 1; //Ground state C-hadron
256  //--------------
257  const xAOD::TruthParticle * parent2=getPreviousParent(parent1);
258  if(!parent2 || !parent2->isHadron()) return 0;
259  if(parent2->isBottomHadron()) return isExcitedHadron(parent2) ? 0 : 1; //Ground state B-hadron
260  if(parent2->isCharmHadron() && !isExcitedHadron(parent2))return 1; //Ground state C-hadron
261  //--------------
262  const xAOD::TruthParticle * parent3=getPreviousParent(parent2);
263  if(!parent3 || !parent3->isHadron()) return 0;
264  if(parent3->isBottomHadron()) return isExcitedHadron(parent3) ? 0 : 1; //Ground state B-hadron
265  if(parent3->isCharmHadron() && !isExcitedHadron(parent3))return 1; //Ground state C-hadron
266  }
267  return 0;
268  }
269  const xAOD::TruthParticle * NewVrtSecInclusiveTool::getPreviousParent(const xAOD::TruthParticle * child) {
270  if( child->hasProdVtx() ){
271  if( child->prodVtx()->nIncomingParticles()==1 ){
272  return *(child->prodVtx()->incomingParticleLinks())[0];
273  }
274  }
275  return nullptr;
276  }
277  bool NewVrtSecInclusiveTool::isExcitedHadron(const xAOD::TruthParticle * tp) {
278  if( !tp->hasProdVtx() )return false;
279  if( !tp->hasDecayVtx() )return false;
280  Amg::Vector3D pvrt(tp->prodVtx()->x(),tp->prodVtx()->y(),tp->prodVtx()->z());
281  Amg::Vector3D dvrt(tp->decayVtx()->x(),tp->decayVtx()->y(),tp->decayVtx()->z());
282  return Amg::distance(pvrt,dvrt) < 0.001*Gaudi::Units::mm ? true : false;
283  }
284 
285  bool NewVrtSecInclusiveTool::isDisplaced(const xAOD::TrackParticle * TP) {
286  if(!TP)return false;
287  static const SG::ConstAccessor<ElementLink<xAOD::TruthParticleContainer>> truthParticleLinkAcc( "truthParticleLink");
288  if( truthParticleLinkAcc.isAvailable( *TP ) ) {
289  const ElementLink<xAOD::TruthParticleContainer>& tplink = truthParticleLinkAcc( *TP );
290  if( !tplink.isValid() ) return false;
291  if((*tplink)->hasProdVtx() && (*tplink)->prodVtx()->perp() > 1.*Gaudi::Units::mm) return true;
292  }
293  return false;
294  }
295 
296  int NewVrtSecInclusiveTool::getG4Inter(const xAOD::TrackParticle* TP ) {
297  static const SG::ConstAccessor< ElementLink< xAOD::TruthParticleContainer> >
298  truthParticleLinkAcc ( "truthParticleLink");
299  if( truthParticleLinkAcc.isAvailable( *TP ) ) {
300  const ElementLink<xAOD::TruthParticleContainer>& tplink =
301  truthParticleLinkAcc( *TP );
302  if( tplink.isValid() && HepMC::is_simulation_particle((*tplink))) return 1;
303  }
304  return 0;
305  }
306  int NewVrtSecInclusiveTool::getMCPileup(const xAOD::TrackParticle* TP ) {
307  static const SG::ConstAccessor< ElementLink< xAOD::TruthParticleContainer> >
308  truthParticleLinkAcc ( "truthParticleLink");
309  if( truthParticleLinkAcc.isAvailable( *TP ) ) {
310  const ElementLink<xAOD::TruthParticleContainer>& tplink =
311  truthParticleLinkAcc( *TP );
312  if( !tplink.isValid() ) return 1;
313  } else { return 1; }
314  return 0;
315  }
316 
317  int NewVrtSecInclusiveTool::getProdVrtBarcode(const xAOD::TrackParticle * TP, float resolLimit) // FIXME barcode-based
318  {
319  int barVrt{HepMC::UNDEFINED_ID};
320  if(!TP)return barVrt;
321  static const SG::ConstAccessor<ElementLink<xAOD::TruthParticleContainer>> truthParticleLinkAcc( "truthParticleLink");
322  if( truthParticleLinkAcc.isAvailable( *TP ) ) {
323  barVrt=-1;
324  const ElementLink<xAOD::TruthParticleContainer>& tplink = truthParticleLinkAcc( *TP );
325  if( !tplink.isValid() ) return barVrt;
326  const xAOD::TruthParticle *tparticle = (*tplink); // truth particle for TrackPartile
327  if(tparticle->hasProdVtx()){ // truth particle has production vertex
328  barVrt = HepMC::uniqueID(tparticle->prodVtx());
329  if(HepMC::is_simulation_vertex(tparticle->prodVtx())) return barVrt; // Geant4 vertex
330  if(tparticle->prodVtx()->nIncomingParticles()!=1) return barVrt; // Safety!
331  const xAOD::TruthParticle *parent = *(tparticle->prodVtx()->incomingParticleLinks())[0];
332  Amg::Vector3D vpos0(tparticle->prodVtx()->x(),tparticle->prodVtx()->y(),tparticle->prodVtx()->z()); //Truth vertex position
333  //
334  //-- Now check parent
335  if(!parent->isHadron()) return barVrt; // Parent is parton not particle! Stop.
336  if(!parent->hasProdVtx()) return barVrt; // Parent particle doesn't have production vertex
337  Amg::Vector3D vpos1(parent->prodVtx()->x(),parent->prodVtx()->y(),parent->prodVtx()->z()); //Truth vertex position
338  if( Amg::distance(vpos0,vpos1) > resolLimit) return barVrt; // Parent vertex is far
339  barVrt = HepMC::uniqueID(parent->prodVtx()); // Else use parent vertex as reference
340  if(HepMC::is_simulation_vertex(parent->prodVtx())) return barVrt; // Geant4 vertex
341  if(parent->prodVtx()->nIncomingParticles()!=1) return barVrt; // Not a decay vertex
342  const xAOD::TruthParticle *grandparent = *(parent->prodVtx()->incomingParticleLinks())[0];
343  //
344  //-- Now check grandparent
345  if(!grandparent->isHadron()) return barVrt; // Grandparent is parton not particle! Stop.
346  if(!grandparent->hasProdVtx()) return barVrt; // Parent particle doesn't have production vertex
347  Amg::Vector3D vpos2(grandparent->prodVtx()->x(),grandparent->prodVtx()->y(),grandparent->prodVtx()->z()); //Truth vertex position
348  if( Amg::distance(vpos0,vpos2) > resolLimit) return barVrt; // Grandparent vertex is far
349  barVrt = HepMC::uniqueID(grandparent->prodVtx()); // Use grandparent vertex as reference
350  if(HepMC::is_simulation_vertex(grandparent->prodVtx())) return barVrt; // Geant4 vertex
351  if(grandparent->prodVtx()->nIncomingParticles()!=1) return barVrt; // Not a decay vertex
352  const xAOD::TruthParticle *biggrandparent = *(grandparent->prodVtx()->incomingParticleLinks())[0];
353  //
354  //-- Now check biggrandparent
355  if(!biggrandparent->isHadron()) return barVrt; // BigGrandparent is parton not particle! Stop.
356  if(!biggrandparent->hasProdVtx()) return barVrt; // Parent particle doesn't have production vertex
357  Amg::Vector3D vpos3(biggrandparent->prodVtx()->x(),biggrandparent->prodVtx()->y(),biggrandparent->prodVtx()->z());
358  if( Amg::distance(vpos0,vpos3) > resolLimit) return barVrt; // Grandparent vertex is far
359  barVrt = HepMC::uniqueID(biggrandparent->prodVtx()); // Use grandparent vertex as reference
360  }
361  }
362  return barVrt;
363  }
364 
365  bool NewVrtSecInclusiveTool::checkTrue2TrVrt(const xAOD::TrackParticle * TP1, const xAOD::TrackParticle * TP2, float nearCut)
366  {
367  if( !TP1 || !TP2 )return false;
368  static const SG::ConstAccessor<ElementLink<xAOD::TruthParticleContainer>> truthParticleLinkAcc( "truthParticleLink");
369  if( !truthParticleLinkAcc.isAvailable(*TP1) ) return false;
370  if( !truthParticleLinkAcc.isAvailable(*TP2) ) return false;
371  const ElementLink<xAOD::TruthParticleContainer>& tplink1 = truthParticleLinkAcc( *TP1 );
372  if( !tplink1.isValid() ) return false;
373  const ElementLink<xAOD::TruthParticleContainer>& tplink2 = truthParticleLinkAcc( *TP2 );
374  if( !tplink2.isValid() ) return false;
375  const xAOD::TruthParticle *tpart1 = (*tplink1); // truth particle for TrackPartile1
376  const xAOD::TruthParticle *tpart2 = (*tplink2); // truth particle for TrackPartile2
377  //-----
378  if(!tpart1->hasProdVtx()) return false;
379  if(!tpart2->hasProdVtx()) return false;
380  if(tpart1->prodVtx()->nOutgoingParticles()<2) return false;
381  if(tpart2->prodVtx()->nOutgoingParticles()<2) return false;
382  //-----
383  Amg::Vector3D vpos1(tpart1->prodVtx()->x(),tpart1->prodVtx()->y(),tpart1->prodVtx()->z());
384  Amg::Vector3D vpos2(tpart2->prodVtx()->x(),tpart2->prodVtx()->y(),tpart2->prodVtx()->z());
385  return Amg::distance(vpos1,vpos2)<nearCut ? true : false;
386  }
387 
388 
389  double NewVrtSecInclusiveTool::distToMatLayerSignificance(Vrt2Tr & Vrt) const
390  {
391  const EventContext& ctx = Gaudi::Hive::currentContext();
392  if(Vrt.fitVertex.perp()<20.) return 1.e9;
393  double normP=1./Vrt.momentum.P();
394  Amg::Vector3D momentumP(Vrt.momentum.Px()*normP,Vrt.momentum.Py()*normP,Vrt.momentum.Pz()*normP);
395  Amg::Vector3D momentumN=-momentumP;
396 
397  const Trk::Layer * someLayer = nullptr;
398  const Trk::Layer * nextLayerP = nullptr;
399  const Trk::Layer * nextLayerN = nullptr;
400  const auto *volume = m_extrapolator->trackingGeometry()->lowestTrackingVolume(Vrt.fitVertex);
401  someLayer = volume->associatedLayer(Vrt.fitVertex);
402  const auto *material = someLayer->layerMaterialProperties();
403  if(material){
404  nextLayerP=someLayer;
405  } else {
406  nextLayerP = someLayer->nextLayer(Vrt.fitVertex,momentumP);
407  if(nextLayerP){ if(!nextLayerP->layerMaterialProperties())nextLayerP=nullptr; }
408  nextLayerN = someLayer->nextLayer(Vrt.fitVertex,momentumN);
409  if(nextLayerN){ if(!nextLayerN->layerMaterialProperties())nextLayerN=nullptr; }
410  }
411  momentumP *= 1.e5; //100GeV to have straight trajectory
412  double charge = 1.;
413  const Trk::Perigee pseudoVrtPart(Vrt.fitVertex, momentumP, charge, Vrt.fitVertex);
414 
415  const Trk::TrackParameters * extrapParP=nullptr; //along momentum
416  const Trk::TrackParameters * extrapParN=nullptr; //backward
417  if(nextLayerP){ extrapParP = m_extrapolator->extrapolate(ctx, pseudoVrtPart,
418  nextLayerP->surfaceRepresentation(), Trk::anyDirection, false, Trk::nonInteractingMuon).release();}
419  if(nextLayerN){ extrapParN = m_extrapolator->extrapolate(ctx, pseudoVrtPart,
420  nextLayerN->surfaceRepresentation(), Trk::anyDirection, false, Trk::nonInteractingMuon).release();}
421 
422  float distanceP=1.e9, distanceN=1.e9;
423  if(extrapParP)distanceP=PntPntDist(extrapParP->position(), Vrt.fitVertex);
424  if(extrapParN)distanceN=PntPntDist(extrapParN->position(), Vrt.fitVertex);
425  if(distanceP==1.e9 && distanceN==1.e9) return 1.e9;
426 
427  //std::pair<const Trk::TrackParameters*,const Trk::Layer*> next=
428  // m_extrapolator->extrapolateToNextActiveLayer(pseudoVrtPart,Trk::anyDirection,true,Trk::pion) ;
429 
430  double signif=1.e9;
431  std::vector<double> pntCovar={1.e-2,0.,1.e-2,0.,0.,4.e-2};
432  if(distanceP<distanceN)signif=vrtVrtDist(Vrt.fitVertex, Vrt.errorMatrix, extrapParP->position(), pntCovar);
433  else signif=vrtVrtDist(Vrt.fitVertex, Vrt.errorMatrix, extrapParN->position(), pntCovar);
434  delete extrapParP;
435  delete extrapParN;
436  return signif;
437  }
438 
439  std::vector<double> NewVrtSecInclusiveTool::estimVrtPos( int nTrk, std::deque<long int> &selTrk, std::map<long int, std::vector<double>> & vrt)
440  {
441  std::vector<double> estimation(3,0.);
442  int ntsel=selTrk.size();
443  for( int i=0; i<ntsel-1; i++){
444  for( int j=i+1; j<ntsel; j++){
445  int k = selTrk[i]<selTrk[j] ? selTrk[i]*nTrk+selTrk[j] : selTrk[j]*nTrk+selTrk[i];
446  estimation[0]+=vrt.at(k)[0];
447  estimation[1]+=vrt[k][1];
448  estimation[2]+=vrt[k][2];
449  } }
450  estimation[0] /= ntsel*(ntsel-1)/2;
451  estimation[1] /= ntsel*(ntsel-1)/2;
452  estimation[2] /= ntsel*(ntsel-1)/2;
453  return estimation;
454  }
455 
456 } //end namespace
Trk::pixelEndCap2
@ pixelEndCap2
Definition: Tracking/TrkEvent/TrkTrackSummary/TrkTrackSummary/TrackSummary.h:239
xAOD::TruthVertex_v1::nOutgoingParticles
size_t nOutgoingParticles() const
Get the number of outgoing particles.
Trk::anyDirection
@ anyDirection
Definition: PropDirection.h:22
Rec::NewVrtSecInclusiveTool::getIdHF
int getIdHF(const xAOD::TrackParticle *TP) const
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:242
Rec::NewVrtSecInclusiveTool::vrtVrtDist2D
static double vrtVrtDist2D(const xAOD::Vertex &primVrt, const Amg::Vector3D &secVrt, const std::vector< double > &vrtErr, double &signif)
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:93
xAOD::Vertex_v1::x
float x() const
Returns the x position.
Rec::NewVrtSecInclusiveTool::getIBLHit
static int getIBLHit(const xAOD::TrackParticle *Part)
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:214
Rec::NewVrtSecInclusiveTool::isDisplaced
static bool isDisplaced(const xAOD::TrackParticle *tp)
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:285
Rec::NewVrtSecInclusiveTool::getBLHit
static int getBLHit(const xAOD::TrackParticle *Part)
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:223
test_pyathena.px
px
Definition: test_pyathena.py:18
python.SystemOfUnits.mm
float mm
Definition: SystemOfUnits.py:98
xAOD::uint8_t
uint8_t
Definition: Muon_v1.cxx:558
xAOD::TruthParticle_v1::isCharmHadron
bool isCharmHadron() const
Determine if the PID is that of a c-hadron.
xAOD::uint32_t
setEventNumber uint32_t
Definition: EventInfo_v1.cxx:127
Trk::ParametersBase::position
const Amg::Vector3D & position() const
Access method for the position.
Trk::ParametersT
Dummy class used to allow special convertors to be called for surfaces owned by a detector element.
Definition: EMErrorDetail.h:25
sincos.h
Helper to simultaneously calculate sin and cos of the same angle.
ParticleTest.tp
tp
Definition: ParticleTest.py:25
AthCommonMsg< AlgTool >::msgLvl
bool msgLvl(const MSG::Level lvl) const
Definition: AthCommonMsg.h:30
Layer.h
Rec::NewVrtSecInclusiveTool::getPixelDiscs
static void getPixelDiscs(const xAOD::TrackParticle *Part, int &d0Hit, int &d1Hit, int &d2Hit)
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:233
xAOD::expectInnermostPixelLayerHit
@ expectInnermostPixelLayerHit
Do we expect a 0th-layer barrel hit for this track?
Definition: TrackingPrimitives.h:237
NeutralParameters.h
SG::ConstAccessor
Helper class to provide constant type-safe access to aux data.
Definition: ConstAccessor.h:55
TrigInDetValidation_Base.test
test
Definition: TrigInDetValidation_Base.py:142
Rec::NewVrtSecInclusiveTool::momAtVrt
ROOT::Math::PxPyPzEVector momAtVrt(const std::vector< double > &inpTrk) const
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:201
AmgSymMatrix
#define AmgSymMatrix(dim)
Definition: EventPrimitives.h:50
xAOD::TruthVertex_v1::y
float y() const
Vertex y displacement.
Rec::NewVrtSecInclusiveTool::Vrt2Tr
Definition: NewVrtSecInclusiveTool.h:298
Rec::NewVrtSecInclusiveTool::massV0
static double massV0(const std::vector< std::vector< double > > &TrkAtVrt, double massP, double massPi)
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:177
TrkVKalVrtFitter.h
Rec::NewVrtSecInclusiveTool::estimVrtPos
static std::vector< double > estimVrtPos(int nTrk, std::deque< long int > &selTrk, std::map< long int, std::vector< double >> &vrt)
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:439
Rec::NewVrtSecInclusiveTool::Vrt2Tr::fitVertex
Amg::Vector3D fitVertex
Definition: NewVrtSecInclusiveTool.h:299
Trk::pixelEndCap0
@ pixelEndCap0
three pixel discs (on each side)
Definition: Tracking/TrkEvent/TrkTrackSummary/TrkTrackSummary/TrackSummary.h:237
python.changerun.kk
list kk
Definition: changerun.py:39
Rec::NewVrtSecInclusiveTool::getMCPileup
static int getMCPileup(const xAOD::TrackParticle *TP)
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:306
xAOD::expectNextToInnermostPixelLayerHit
@ expectNextToInnermostPixelLayerHit
Do we expect a 1st-layer barrel hit for this track?
Definition: TrackingPrimitives.h:248
Trk::Layer::surfaceRepresentation
virtual const Surface & surfaceRepresentation() const =0
Transforms the layer into a Surface representation for extrapolation.
CxxUtils::sincos::cs
double cs
Definition: sincos.h:54
HepMC::is_simulation_particle
bool is_simulation_particle(const T &p)
Method to establish if a particle (or barcode) was created during the simulation (TODO update to be s...
Definition: MagicNumbers.h:354
Rec::NewVrtSecInclusiveTool::checkTrue2TrVrt
static bool checkTrue2TrVrt(const xAOD::TrackParticle *TP1, const xAOD::TrackParticle *TP2, float nearCut=0.1)
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:365
lumiFormat.i
int i
Definition: lumiFormat.py:85
Rec
Name: MuonSpContainer.h Package : offline/Reconstruction/MuonIdentification/muonEvent.
Definition: FakeTrackBuilder.h:10
Trk::pixelEndCap1
@ pixelEndCap1
Definition: Tracking/TrkEvent/TrkTrackSummary/TrkTrackSummary/TrackSummary.h:238
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
xAOD::TruthParticle_v1
Class describing a truth particle in the MC record.
Definition: TruthParticle_v1.h:37
Rec::NewVrtSecInclusiveTool::m_massPi
double m_massPi
Definition: NewVrtSecInclusiveTool.h:195
xAOD::TruthParticle_v1::hasProdVtx
bool hasProdVtx() const
Check for a production vertex on this particle.
Definition: TruthParticle_v1.cxx:69
HepMC::uniqueID
int uniqueID(const T &p)
Definition: MagicNumbers.h:116
Amg::pz
@ pz
Definition: GeoPrimitives.h:40
Rec::NewVrtSecInclusiveTool::getProdVrtBarcode
static int getProdVrtBarcode(const xAOD::TrackParticle *TP, float resolLimit=0.1)
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:317
test_pyathena.parent
parent
Definition: test_pyathena.py:15
Trk::Layer::nextLayer
const Layer * nextLayer(const Amg::Vector3D &gp, const Amg::Vector3D &udir) const
getting the next/previous Layer if registered - unit for direction vector required
Definition: Layer.cxx:161
HepMC::is_simulation_vertex
bool is_simulation_vertex(const T &v)
Method to establish if the vertex was created during simulation (TODO migrate to be based on status).
Definition: MagicNumbers.h:360
TrackSummary.h
Trk::ParametersBase
Definition: ParametersBase.h:55
xAOD::Vertex_v1::z
float z() const
Returns the z position.
HepMC::UNDEFINED_ID
constexpr int UNDEFINED_ID
Definition: MagicNumbers.h:56
xAOD::TruthParticle_v1::isBottomHadron
bool isBottomHadron() const
Determine if the PID is that of a b-hadron.
xAOD::TruthParticle_v1::prodVtx
const TruthVertex_v1 * prodVtx() const
The production vertex of this particle.
Definition: TruthParticle_v1.cxx:75
Amg::py
@ py
Definition: GeoPrimitives.h:39
xAOD::numberOfNextToInnermostPixelLayerHits
@ numberOfNextToInnermostPixelLayerHits
these are the hits in the 1st pixel barrel layer
Definition: TrackingPrimitives.h:249
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
MagicNumbers.h
Rec::NewVrtSecInclusiveTool::getPreviousParent
static const xAOD::TruthParticle * getPreviousParent(const xAOD::TruthParticle *child)
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:269
charge
double charge(const T &p)
Definition: AtlasPID.h:986
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
CxxUtils::sincos::sn
double sn
Definition: sincos.h:51
NewVrtSecInclusiveTool.h
xAOD::TruthVertex_v1::x
float x() const
Vertex x displacement.
makeTRTBarrelCans.dy
tuple dy
Definition: makeTRTBarrelCans.py:21
TrackingVolume.h
Rec::NewVrtSecInclusiveTool::Vrt2Tr::errorMatrix
std::vector< double > errorMatrix
Definition: NewVrtSecInclusiveTool.h:302
Rec::NewVrtSecInclusiveTool::m_extrapolator
ToolHandle< Trk::IExtrapolator > m_extrapolator
Definition: NewVrtSecInclusiveTool.h:184
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
xAOD::Vertex_v1
Class describing a Vertex.
Definition: Vertex_v1.h:42
Rec::NewVrtSecInclusiveTool::vrtVrtDist
static double vrtVrtDist(const xAOD::Vertex &primVrt, const Amg::Vector3D &secVrt, const std::vector< double > &vrtErr, double &signif)
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:60
python.Constants.INFO
int INFO
Definition: Control/AthenaCommon/python/Constants.py:15
xAOD::TruthVertex_v1::nIncomingParticles
size_t nIncomingParticles() const
Get the number of incoming particles.
GeoPrimitivesHelpers.h
xAOD::TruthVertex_v1::z
float z() const
Vertex longitudinal distance along the beam line form the origin.
makeTRTBarrelCans.dx
tuple dx
Definition: makeTRTBarrelCans.py:20
xAOD::Vertex_v1::y
float y() const
Returns the y position.
AthCommonMsg< AlgTool >::msg
MsgStream & msg() const
Definition: AthCommonMsg.h:24
DetType::Part
Part
Definition: DetType.h:14
Trk::Layer::layerMaterialProperties
const LayerMaterialProperties * layerMaterialProperties() const
getting the LayerMaterialProperties including full/pre/post update
Rec::NewVrtSecInclusiveTool::getG4Inter
static int getG4Inter(const xAOD::TrackParticle *TP)
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:296
CxxUtils::sincos
Helper to simultaneously calculate sin and cos of the same angle.
Definition: sincos.h:39
Rec::NewVrtSecInclusiveTool::distToMatLayerSignificance
double distToMatLayerSignificance(Vrt2Tr &Vrt) const
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:389
SG::ConstAccessor::isAvailable
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
TrackingGeometry.h
Rec::NewVrtSecInclusiveTool::projSV_PV
static double projSV_PV(const Amg::Vector3D &SV, const xAOD::Vertex &PV, const TLorentzVector &Direction)
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:54
xAOD::TrackParticle_v1
Class describing a TrackParticle.
Definition: TrackParticle_v1.h:43
Trk::nonInteractingMuon
@ nonInteractingMuon
Definition: ParticleHypothesis.h:39
ConstAccessor.h
Helper class to provide constant type-safe access to aux data.
Rec::NewVrtSecInclusiveTool::vrtRadiusError
static double vrtRadiusError(const Amg::Vector3D &secVrt, const std::vector< double > &vrtErr)
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:160
Rec::NewVrtSecInclusiveTool::printWrkSet
void printWrkSet(const std::vector< WrkVrt > *WrkSet, const std::string &name) const
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:27
Rec::NewVrtSecInclusiveTool::isExcitedHadron
static bool isExcitedHadron(const xAOD::TruthParticle *tp)
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:277
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
Rec::NewVrtSecInclusiveTool::Vrt2Tr::momentum
TLorentzVector momentum
Definition: NewVrtSecInclusiveTool.h:300
fitman.k
k
Definition: fitman.py:528
xAOD::numberOfInnermostPixelLayerHits
@ numberOfInnermostPixelLayerHits
these are the hits in the 0th pixel barrel layer
Definition: TrackingPrimitives.h:238
xAOD::TruthParticle_v1::isHadron
bool isHadron() const
Whether the particle is a hadron.
Trk::Layer
Definition: Layer.h:72
xAOD::TruthVertex_v1::incomingParticleLinks
const TPLinks_t & incomingParticleLinks() const
Get all the incoming particles.
Rec::NewVrtSecInclusiveTool::PntPntDist
static double PntPntDist(const Amg::Vector3D &Vrt1, const Amg::Vector3D &Vrt2)
Definition: Reconstruction/VKalVrt/NewVrtSecInclusiveTool/src/Utilities.cxx:149
Generated on Thu Sep 4 2025 21:23:22 for ATLAS Offline Software by   doxygen 1.8.18