ATLAS Offline Software
VtDeriv.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/VtDeriv.h"
6 #include "TrkVKalVrtCore/CommonPars.h"
7 #include "TrkVKalVrtCore/Propagator.h"
8 #include "TrkVKalVrtCore/Matrix.h"
9 #include "TrkVKalVrtCore/Utilities.h"
10 #include "TrkVKalVrtCore/TrkVKalVrtCore.h"
11 #include "TrkVKalVrtCore/VKalVrtBMag.h"
12 #include <cmath>
13 
14 namespace Trk {
15 
16 void vpderiv(bool UseTrackErr, long int Charge, const double *pari0, double *covi, double *vrtref, double *covvrtref,
17  double *drdpar, double *dwgt, double *rv0, VKalVrtControl * FitCONTROL)
18 {
19  /* Initialized data */
20 
21  double dshft[6] = { 1.e-4, 1.e-4, 1.e-4, 1.e-4, 1.e-4, 1.e-4}; //step for each parameter
22  /* Symmetric version for derivatives calculation*/
23  double DrvCoef[4] = {-2., -1., 1.,2.};
24 
25  /* Local variables */
26  double pari[6], covd[15], dcov[3], rvec[50]; /* was [2][6*4+1] */
27  double paro[5];
28  int j, ij, ip, ipp, id=0;
29  double dwgt0[3]={0.,0.,0.};
30  //double deriv[6], dchi2[4*6+1]; //VK for debugging
31  double cs, pp, sn, pt, rho;
32  double covdtr[15], ctg, par[5], cnv[36]; /* was [6][6] */
33 
34 #define rvec_ref(a_1,a_2) rvec[(a_2)*2 + (a_1) - 1]
35 #define drdpar_ref(a_1,a_2) drdpar[(a_2)*2 + (a_1)]
36 #define cnv_ref(a_1,a_2) cnv[(a_2)*6 + (a_1) - 7]
37 
38 /* ---------------------------------------------------------- */
39 /* Subroutine for derivative calculations */
40 /* Input: */
41 /* Charge - charge of track ( +1,0,-1 ) */
42 /* PARI0(1:3) - Vertex of particle */
43 /* PARI0(4:6) - Momentum of particle */
44 /* COVI(21) - symmetric covariance matrix */
45 /* BMAG - magnetic field */
46 /* UseTrackErr =false Fitted track is assumed errorless */
47 /* only vertex error is used */
48 /* may be used for constraint */
49 /* =true Fitted track error is added */
50 /* Output: */
51 /* Deriv(6) - Chi2 derivatives for debugging */
52 /* drdpar - derivatives */
53 /* rv0[2] - track impact vs vertex */
54 /* dwgt[3] - track impact covarinance */
55 /* ---------------------------------------------------------- */
56  /* Parameter adjustments */
57  drdpar -= 3;
58  --covvrtref;
59 
60  /* Function Body */
61 /* --------------------- */
62  double constB =FitCONTROL->vk_forcft.localbmag * vkalMagCnvCst ;
63 
64  for (ip = 0; ip <= 4*6; ++ip) { // Number of points * Number of parameters
65 /* -- Input parameters */
66  pari[0] = pari0[0];
67  pari[1] = pari0[1];
68  pari[2] = pari0[2];
69  pari[3] = pari0[3];
70  pari[4] = pari0[4];
71  pari[5] = pari0[5];
72 // ipp - parameter number, id - number of step
73  ipp = ip; while(ipp>6)ipp-=6;
74  if(ip>=1 && ip<=6 ) id=0;
75  if(ip>=7 && ip<=12) id=1;
76  if(ip>=13 && ip<=18) id=2;
77  if(ip>=19 && ip<=24) id=3;
78  if (ip != 0) pari[ipp - 1] *= (dshft[ipp-1]*DrvCoef[id] + 1.);
79 /* -- */
80  pt = sqrt(pari[3]*pari[3] + pari[4]*pari[4]);
81  pp = pt*pt + pari[5]*pari[5];
82  cs = pari[3] / pt;
83  sn = pari[4] / pt;
84  ctg = pari[5] / pt;
85  rho = Charge * constB / pt; // Sign of charge coinside with sign of curvature!!!
86 /* -- Output parameters */
87  par[0] = 0.;
88  par[1] = 0.;
89  par[2] = acos(pari[5] / sqrt(pp));
90  if(par[2]<1.e-5) par[2]=1.e-5;
91  if(par[2]>M_PI-1.e-5) par[2]=M_PI-1.e-5;
92  par[3] = atan2(pari[4], pari[3]);
93  par[4] = rho;
94  if (Charge == 0) par[4] = constB / pt;
95 /* -- */
96  cnv_ref(1, 1) = sn;
97  cnv_ref(2, 1) = -cs;
98  cnv_ref(3, 1) = 0.;
99  cnv_ref(4, 1) = 0.;
100  cnv_ref(5, 1) = 0.;
101  cnv_ref(6, 1) = 0.;
102 /* -- */
103  cnv_ref(1, 2) = -cs * ctg;
104  cnv_ref(2, 2) = -sn * ctg;
105  cnv_ref(3, 2) = 1.;
106  cnv_ref(4, 2) = 0.;
107  cnv_ref(5, 2) = 0.;
108  cnv_ref(6, 2) = 0.;
109 /* -- */
110  cnv_ref(1, 3) = 0.;
111  cnv_ref(2, 3) = 0.;
112  cnv_ref(3, 3) = 0.;
113  cnv_ref(4, 3) = pari[5] * pari[3] / pp / pt;
114  cnv_ref(5, 3) = pari[5] * pari[4] / pp / pt;
115  cnv_ref(6, 3) = -pt / pp;
116 /* -- */
117  cnv_ref(1, 4) = 0.;
118  cnv_ref(2, 4) = 0.;
119  if (Charge)cnv_ref(1, 4) = -cs * rho; // For charged tracks only
120  if (Charge)cnv_ref(2, 4) = -sn * rho; //
121  cnv_ref(3, 4) = 0.;
122  cnv_ref(4, 4) = -pari[4] / (pt * pt);
123  cnv_ref(5, 4) = pari[3] / (pt * pt);
124  cnv_ref(6, 4) = 0.;
125 /* -- */
126  cnv_ref(1, 5) = 0.;
127  cnv_ref(2, 5) = 0.;
128  cnv_ref(3, 5) = 0.;
129  cnv_ref(4, 5) = -pari[3] * rho / (pt * pt);
130  cnv_ref(5, 5) = -pari[4] * rho / (pt * pt);
131  if ( Charge == 0) {cnv_ref(4, 5) = -pari[3] * constB / (pt*pt*pt);}
132  if ( Charge == 0) {cnv_ref(5, 5) = -pari[4] * constB / (pt*pt*pt);}
133  cnv_ref(6, 5) = 0.;
134 /* -- not needed really but... */
135  cnv_ref(1, 6) = cnv_ref(2, 6) = cnv_ref(3, 6) = cnv_ref(4, 6) = cnv_ref(5, 6) = cnv_ref(6, 6) = 0.;
136 
137  tdasatVK(cnv, &covi[0], covd, 5, 6);
138 
139 /* -- Translation to New Reference Point */
140  Trk::vkalPropagator::Propagate(-999, Charge, par, covd, pari, vrtref, paro, covdtr, FitCONTROL);
141 
142 
143 
144 /* -- Now impact significance */
145 /* X=paro(1)*sn, Y=-paro(1)*cs */
146  sn = sin(paro[3]);
147  cs = cos(paro[3]);
148 /* -- New error version */
149  cnv_ref(1, 1) = -sn;
150  cnv_ref(1, 2) = cs;
151  cnv_ref(1, 3) = 0.;
152  cnv_ref(2, 1) = sn / tan(paro[2]);
153  cnv_ref(2, 2) = cs / tan(paro[2]);
154  cnv_ref(2, 3) = -1.;
155  dcov[0] = 0.;
156  dcov[1] = 0.;
157  dcov[2] = 0.;
158  for (int i = 1; i <= 3; ++i) {
159  for (j = 1; j <= 3; ++j) {
160  ij = i > j ? i : j;
161  ij = ij * (ij - 1) / 2 + ( i < j ? i : j );
162  dcov[0] += cnv_ref(1, i) * cnv_ref(1, j) * covvrtref[ij];
163  dcov[2] += cnv_ref(2, i) * cnv_ref(2, j) * covvrtref[ij];
164  dcov[1] += cnv_ref(1, i) * cnv_ref(2, j) * covvrtref[ij];
165  }
166  }
167 /* --------------------------------------------------------------- */
168  if( UseTrackErr ){ // Add track error if needed
169  dcov[0] += covdtr[0];
170  dcov[1] += covdtr[1];
171  dcov[2] += covdtr[2];
172  }
173 /*---------------------------------------------------------------- */
174 // Weight matrix and chi2 for given shift
175 //
176  int jerr=cfdinv(dcov, &dwgt[0], -2); if(jerr){jerr=cfdinv(dcov, &dwgt[0], 2); if(jerr){dwgt[0]=dwgt[2]=1.e6; dwgt[1]=0.;}};
177  //dchi2[ip] = sqrt(fabs(dwgt[0]*paro[0]*paro[0] + 2.*dwgt[1]*paro[0]*paro[1] + dwgt[2]*paro[1]*paro[1]));
178  rvec_ref(1, ip) = paro[0];
179  rvec_ref(2, ip) = paro[1];
180  if (ip == 0) { //Save weight matrix for output
181  dwgt0[0] = dwgt[0];
182  dwgt0[1] = dwgt[1];
183  dwgt0[2] = dwgt[2];
184  }
185  }
186 /* ---------------------------------------------------------------- */
187 // Now calculate derivatives themselves
188 //
189  for (int i = 1; i <= 6; ++i) {
190  //deriv[i-1] = (dchi2[i] -8.*dchi2[i+6] +8.*dchi2[i+12]-dchi2[i+18] ) / (12.*dshft[i-1]*pari0[i-1]);
191  drdpar_ref(1,i) = (rvec_ref(1,i) -8.*rvec_ref(1,i+6) +8.*rvec_ref(1,i+12) -rvec_ref(1,i+18)) / (12.*dshft[i-1]*pari0[i-1]);
192  drdpar_ref(2,i) = (rvec_ref(2,i) -8.*rvec_ref(2,i+6) +8.*rvec_ref(2,i+12) -rvec_ref(2,i+18)) / (12.*dshft[i-1]*pari0[i-1]);
193  }
194 //std::cout<<deriv[0]<<", "<<deriv[1]<<", "<<deriv[2]<<", "<<deriv[3]<<", "<<deriv[4]<<", "<<deriv[5]<<'\n';
195 //std::cout<<drdpar_ref(1, 1)<<", "<<drdpar_ref(1, 2)<<", "<<drdpar_ref(1, 3)<<", "<<drdpar_ref(1, 4)<<
196 // ", "<<drdpar_ref(1, 5)<<", "<<drdpar_ref(1, 6)<<'\n';
197 //std::cout<<drdpar_ref(2, 1)<<", "<<drdpar_ref(2, 2)<<", "<<drdpar_ref(2, 3)<<", "<<drdpar_ref(1, 4)<<
198 // ", "<<drdpar_ref(2, 5)<<", "<<drdpar_ref(2, 6)<<'\n';
199 //std::cout<<"---------------------"<<'\n';
200 //
201 // Output parameters
202 //
203 
204  rv0[0] = rvec_ref(1, 0);
205  rv0[1] = rvec_ref(2, 0);
206  dwgt[0] = dwgt0[0];
207  dwgt[1] = dwgt0[1];
208  dwgt[2] = dwgt0[2];
209 
210 //std::cout<<" VpDerivR,Z="<<rv0[0]<<", "<<rv0[1]<<'\n';
211 //std::cout<<" Ref="<<vrtref[0]<<", "<<vrtref[1]<<", "<<vrtref[2]<<'\n';
212 #undef cnv_ref
213 #undef drdpar_ref
214 #undef rvec_ref
215 }
216 } /* End of namespace */
Propagator.h
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
CommonPars.h
test_pyathena.pt
pt
Definition: test_pyathena.py:11
M_PI
#define M_PI
Definition: ActiveFraction.h:11
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
VKalVrtBMag.h
VtDeriv.h
Trk::cfdinv
int cfdinv(double *cov, double *wgt, long int NI)
Definition: Tracking/TrkVertexFitter/TrkVKalVrtCore/src/Matrix.cxx:497
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
lumiFormat.i
int i
Definition: lumiFormat.py:85
vkalMagCnvCst
#define vkalMagCnvCst
Definition: CommonPars.h:23
cnv_ref
#define cnv_ref(a_1, a_2)
TrkVKalVrtCore.h
Matrix.h
find_tgc_unfilled_channelids.ip
ip
Definition: find_tgc_unfilled_channelids.py:3
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
Trk::tdasatVK
void tdasatVK(const double *Der, const double *CovI, double *CovF, long int M, long int N)
Definition: Tracking/TrkVertexFitter/TrkVKalVrtCore/src/Utilities.cxx:218
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
rvec_ref
#define rvec_ref(a_1, a_2)
id
SG::auxid_t id
Definition: Control/AthContainers/Root/debug.cxx:227
createCoolChannelIdFile.par
par
Definition: createCoolChannelIdFile.py:29
Trk::ForCFT::localbmag
double localbmag
Definition: ForCFT.h:28
drdpar_ref
#define drdpar_ref(a_1, a_2)
Trk::VKalVrtControl
Definition: TrkVKalVrtCore.h:44
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
Utilities.h
Trk::VKalVrtControl::vk_forcft
ForCFT vk_forcft
Definition: TrkVKalVrtCore.h:91
fitman.rho
rho
Definition: fitman.py:532
Trk::vpderiv
void vpderiv(bool UseTrackErr, long int Charge, const double *pari0, double *covi, double *vrtref, double *covvrtref, double *drdpar, double *dwgt, double *rv0, VKalVrtControl *FitCONTROL)
Definition: VtDeriv.cxx:16
Generated on Mon Dec 23 2024 21:22:58 for ATLAS Offline Software by   doxygen 1.8.18