ATLAS Offline Software
cfImp.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "TrkVKalVrtCore/cfImp.h"
6 #include "TrkVKalVrtCore/Matrix.h"
7 #include "TrkVKalVrtCore/Propagator.h"
8 #include "TrkVKalVrtCore/TrkVKalVrtCore.h"
9 #include <cmath>
10 
11 namespace {
12 //anonymous namespace.
13 //For internal methods
14 void cfClstPnt(double *par, const double *Vrt, double *ClstPnt) noexcept {
15  double e[3]; // Track direction at perigee
16  e[0] = cos(par[3]);
17  e[1] = sin(par[3]);
18  e[2] = 1. / tan(par[2]);
19  double e2 = 1. + e[2] * e[2]; // vector length squared
20 
21  double Per[3]; // Perigee position
22  Per[0] = sin(par[3]) * par[0];
23  Per[1] = -cos(par[3]) * par[0];
24  Per[2] = par[1];
25 
26  double u = (Vrt[0] - Per[0]) * e[0] + (Vrt[1] - Per[1]) * e[1] +
27  (Vrt[2] - Per[2]) * e[2];
28  u = u / e2;
29 
30  ClstPnt[0] = Per[0] + u * e[0];
31  ClstPnt[1] = Per[1] + u * e[1];
32  ClstPnt[2] = Per[2] + u * e[2];
33 }
34 }
35 
36 
37 namespace Trk {
38 
39 
40 #define min(a,b) ((a) <= (b) ? (a) : (b))
41 #define max(a,b) ((a) >= (b) ? (a) : (b))
42 
43 void cfimp(long int TrkID, long int ich, int IFL, double *par,
44  const double *err, double *vrt, double *vcov, double *rimp,
45  double *rcov, double *sign , VKalVrtControlBase * FitCONTROL )
46 {
47  double dcov[3], errd[15], paro[5];
48  double dwgt[3], errn[15];
49  int i__, j, ij;
50 
51  double cnv[6] /* was [2][3] */;
52 /* --------------------------------------------------------- */
53 /* Impact parameter calculation */
54 /* Input: */
55 /* ICH -charge(-1,0,1) */
56 /* IFL = 1 contribution from VRT is added */
57 /* IFL =-1 contribution from VRT is subtructed */
58 /* PAR(5),ERR(15) - peregee parameters and error matrix */
59 /* VRT(3),VCOV(6) - vertex and its error matrix */
60 /* */
61 /* SIGNIFICANCE IS CALCULATED FOR PERIGEE POINT BY DEF. */
62 /* (NOT FOR THE CLOSEST POINT!!!) */
63 /* */
64 /* Output: */
65 /* RIMP(1) - impact in XY */
66 /* RIMP(2) - impact in Z */
67 /* RIMP(3) - Theta at new vertex */
68 /* RIMP(4) - Phi at new vertex */
69 /* RIMP(5) - curvature at vertex */
70 /* RCOV(3) - error matrix for RIMP */
71 /* SIGN - impact significance */
72 /* Author: V.Kostyukhin */
73 /* --------------------------------------------------------- */
74  /* Parameter adjustments */
75  --vcov;
76 
77  /* Function Body */
78  for (int ii = 0; ii < 15; ++ii) {errd[ii] = err[ii];}
79 
80 
81  double Ref0[3]={0.,0.,0.}; //base reference point for standard perigee
82 
83  Trk::vkalPropagator::Propagate(TrkID, ich, par, errd, Ref0, vrt, paro, errn, FitCONTROL);
84 
85 //std::cout <<" CFImp new par R,Z="<<paro[0]<<", "<<paro[1]<<'\n';
86 /* ---------- */
87  rimp[0] = paro[0];
88  rimp[1] = paro[1];
89  rimp[2] = paro[2];
90  rimp[3] = paro[3];
91  rimp[4] = paro[4];
92 /* X=paro(1)*sn, Y=-paro(1)*cs */
93  double sn = sin(paro[3]);
94  double cs = cos(paro[3]);
95 /* -- New error version */
96  cnv[0] = -sn;
97  cnv[2] = cs;
98  cnv[4] = 0.;
99  cnv[1] = sn / tan(paro[2]);
100  cnv[3] = cs / tan(paro[2]);
101  cnv[5] = -1.;
102  dcov[0] = 0.;
103  dcov[1] = 0.;
104  dcov[2] = 0.;
105  for (i__ = 1; i__ <= 3; ++i__) {
106  for (j = 1; j <= 3; ++j) {
107  ij = max(i__, j);
108  ij = ij * (ij - 1) / 2 + min(i__, j);
109  dcov[0] += cnv[(i__ << 1) - 2] * cnv[(j << 1) - 2] * vcov[ij];
110  dcov[2] += cnv[(i__ << 1) - 1] * cnv[(j << 1) - 1] * vcov[ij];
111  dcov[1] += cnv[(i__ << 1) - 2] * cnv[(j << 1) - 1] * vcov[ij];
112  }
113  }
114  /* --------------------------------------------------------------- */
115  if (IFL == 1) {
116  rcov[0] = errn[0] + dcov[0];
117  rcov[1] = errn[1] + dcov[1];
118  rcov[2] = errn[2] + dcov[2];
119  } else if (IFL == -1) {
120  rcov[0] = errn[0] - dcov[0];
121  rcov[1] = errn[1] - dcov[1];
122  rcov[2] = errn[2] - dcov[2];
123  } else {
124  rcov[0] = errn[0];
125  rcov[1] = errn[1];
126  rcov[2] = errn[2];
127  }
128  int jerr = cfdinv(rcov, dwgt, -2);
129  if (jerr) {jerr=cfdinv(rcov, dwgt, 2); if(jerr){dwgt[0]=dwgt[2]=1.e6; dwgt[1]=0.;}}
130  (*sign) = sqrt(std::abs(dwgt[0] * rimp[0] * rimp[0] + dwgt[1] * 2. * rimp[0] * rimp[1] +
131  dwgt[2] * rimp[1] * rimp[1]));
132 }
133 
134 
135  void cfimpc(long int TrkID, long int ich, int IFL, double *par,
136  const double *err, double *vrt, double *vcov, double *rimp,
137  double *rcov, double *sign, VKalVrtControlBase * FitCONTROL )
138 {
139  double dcov[3], errd[15], paro[5];
140  double dwgt[3], errn[15], cnv[6]; /* was [2][3] */
141  int i__, j, ij;
142 
143 
144  double cs, sn;
145 /* --------------------------------------------------------- */
146 /* SIGNIFICANCE IS CALCULATED FOR THE CLOSEST POINT NOW!!!*/
147 /* Author: V.Kostyukhin */
148 /* --------------------------------------------------------- */
149  /* Parameter adjustments */
150  --vcov;
151  for (int ii = 0; ii < 15; ++ii) {errd[ii] = err[ii];}
152 
153  double Ref0[3]={0.,0.,0.}; //base reference point for standard perigee
154  Trk::vkalPropagator::Propagate(TrkID, ich, par, errd, Ref0, vrt, paro, errn, FitCONTROL);
155 
156  double tmpVrt[3]={0.,0.,0.}; double ClosestPnt[3];
157  cfClstPnt( paro, tmpVrt, ClosestPnt);
158  paro[0]=sqrt(ClosestPnt[0]*ClosestPnt[0] + ClosestPnt[1]*ClosestPnt[1]);
159  paro[1]=ClosestPnt[2];
160 /* ---------- */
161  rimp[0] = paro[0];
162  rimp[1] = paro[1];
163  rimp[2] = paro[2];
164  rimp[3] = paro[3];
165  rimp[4] = paro[4];
166  sn = sin(paro[3]);
167  cs = cos(paro[3]);
168 /* -- New error version */
169  cnv[0] = -sn;
170  cnv[2] = cs;
171  cnv[4] = 0.;
172  cnv[1] = sn / tan(paro[2]);
173  cnv[3] = cs / tan(paro[2]);
174  cnv[5] = -1.;
175  dcov[0] = 0.;
176  dcov[1] = 0.;
177  dcov[2] = 0.;
178  for (i__ = 1; i__ <= 3; ++i__) {
179  for (j = 1; j <= 3; ++j) {
180  ij = max(i__, j);
181  ij = ij * (ij - 1) / 2 + min(i__, j);
182  dcov[0] += cnv[(i__ << 1) - 2] * cnv[(j << 1) - 2] * vcov[ij];
183  dcov[2] += cnv[(i__ << 1) - 1] * cnv[(j << 1) - 1] * vcov[ij];
184  dcov[1] += cnv[(i__ << 1) - 2] * cnv[(j << 1) - 1] * vcov[ij];
185  }
186  }
187  /* --------------------------------------------------------------- */
188  if (IFL == 1) {
189  rcov[0] = errn[0] + dcov[0];
190  rcov[1] = errn[1] + dcov[1];
191  rcov[2] = errn[2] + dcov[2];
192  } else if (IFL == -1) {
193  rcov[0] = errn[0] - dcov[0];
194  rcov[1] = errn[1] - dcov[1];
195  rcov[2] = errn[2] - dcov[2];
196  } else {
197  rcov[0] = errn[0];
198  rcov[1] = errn[1];
199  rcov[2] = errn[2];
200  }
201  int jerr = cfdinv(rcov, dwgt, -2);
202  if (jerr) {
203  jerr = cfdinv(rcov, dwgt, 2);
204  if (jerr) {
205  dwgt[0] = dwgt[2] = 1.e6;
206  dwgt[1] = 0.;
207  }
208  }
209  (*sign) = sqrt(std::abs(dwgt[0] * rimp[0] * rimp[0] +
210  dwgt[1] * 2. * rimp[0] * rimp[1] +
211  dwgt[2] * rimp[1] * rimp[1]));
212  }
213 
214 } /* end of namespace */
215 
Propagator.h
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
max
#define max(a, b)
Definition: cfImp.cxx:41
Trk::cfimpc
void cfimpc(long int TrkID, long int ich, int IFL, double *par, const double *err, double *vrt, double *vcov, double *rimp, double *rcov, double *sign, VKalVrtControlBase *FitCONTROL)
Definition: cfImp.cxx:135
Trk::cfimp
void cfimp(long int TrkID, long int ich, int IFL, double *par, const double *err, double *vrt, double *vcov, double *rimp, double *rcov, double *sign, VKalVrtControlBase *FitCONTROL)
Definition: cfImp.cxx:43
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
Trk::cfdinv
int cfdinv(double *cov, double *wgt, long int NI)
Definition: Tracking/TrkVertexFitter/TrkVKalVrtCore/src/Matrix.cxx:497
Trk::u
@ u
Enums for curvilinear frames.
Definition: ParamDefs.h:77
Trk::vkalPropagator::Propagate
static void Propagate(long int TrkID, long int Charge, double *ParOld, double *CovOld, double *RefStart, double *RefEnd, double *ParNew, double *CovNew, VKalVrtControlBase *FitControl=0)
Definition: Propagator.cxx:127
dqt_zlumi_pandas.err
err
Definition: dqt_zlumi_pandas.py:182
TrkVKalVrtCore.h
Matrix.h
sign
int sign(int a)
Definition: TRT_StrawNeighbourSvc.h:107
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
min
#define min(a, b)
Definition: cfImp.cxx:40
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
createCoolChannelIdFile.par
par
Definition: createCoolChannelIdFile.py:29
egammaEnergyPositionAllSamples::e2
double e2(const xAOD::CaloCluster &cluster)
return the uncorrected cluster energy in 2nd sampling
Trk::VKalVrtControlBase
Definition: TrkVKalVrtCore.h:26
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
cfImp.h
Generated on Thu Nov 7 2024 21:11:51 for ATLAS Offline Software by   doxygen 1.8.18