134{
135 const double C=0.02999975/1000.0;
136 const double minStep=30.0;
137
138 double J[5][5],Rf[5],AG[5][5],Gf[5][5],
A[5][5];
140
141 bool samePlane=false;
142
143 if(pSB!=nullptr)
144 {
148 samePlane=true;
150 }
151 }
152
153 if(!samePlane) {
154
155 double gP[3],gPi[3],lP[3],gV[3],
a,
b,
c,
s,J0[7][5],
descr,CQ,Ac,Av,Cc;
156 double V[3],
P[3],M[3][3],D[4],Jm[7][7],
157 J1[5][7],gB[3],gBi[3],gBf[3],dBds[3],Buf[5][7],DVx,DVy,DVz;
158 int nStep,nStepMax;
160
161 double sint,
cost,sinf,cosf;
165
166 memset(&J0[0][0],0,sizeof(J0));
167
168 if(pSB!=nullptr)
169 {
174
175 J0[0][0]=
L[0][0];J0[0][1]=
L[0][1];
176 J0[1][0]=
L[1][0];J0[1][1]=
L[1][1];
177 J0[2][0]=
L[2][0];J0[2][1]=
L[2][1];
178 J0[3][2]=-sinf*sint;J0[3][3]=cosf*
cost;
179 J0[4][2]= cosf*sint;J0[4][3]=sinf*
cost;
180 J0[5][3]=-sint;
181 J0[6][4]=1.0;
182 }
183 else
184 {
190 J0[2][1]=1.0;
191 J0[3][2]=-sinf*sint;J0[3][3]=cosf*
cost;
192 J0[4][2]= cosf*sint;J0[4][3]=sinf*
cost;
193 J0[5][3]=-sint;
194 J0[6][4]=1.0;
195 }
196 for(i=0;
i<4;
i++) D[i]=pSE->
getPar(i);
197 for(i=0;
i<3;
i++) gPi[i]=gP[i];
198
200
201 for(i=0;
i<3;
i++) gBi[i]=gB[i];
202
203 c=D[0]*gP[0]+D[1]*gP[1]+D[2]*gP[2]+D[3];
204 b=D[0]*gV[0]+D[1]*gV[1]+D[2]*gV[2];
205 a=0.5*CQ*(gB[0]*(D[1]*gV[2]-D[2]*gV[1])+
206 gB[1]*(D[2]*gV[0]-D[0]*gV[2])+
207 gB[2]*(D[0]*gV[1]-D[1]*gV[0]));
208
210
211 if(descr<0.0)
212 {
213
214 return nullptr;
215 }
216
217 bool useExpansion=true;
219
220 if(fabs(ratio)>0.1)
221 useExpansion = false;
222
223 if(useExpansion) {
226 }
227 else {
228 int signb = (
b<0.0)?-1:1;
229 sl = (-
b+signb*sqrt(descr))/(2*
a);
230 }
231
232 if(fabs(sl)<minStep) nStepMax=1;
233 else
234 {
235 nStepMax=(
int)(fabs(sl)/minStep)+1;
236 }
237 if((nStepMax<0)||(nStepMax>1000))
238 {
239 return nullptr;
240 }
241 Av=sl*CQ;
242 Ac=0.5*sl*Av;
243 DVx=gV[1]*gB[2]-gV[2]*gB[1];
244 DVy=gV[2]*gB[0]-gV[0]*gB[2];
245 DVz=gV[0]*gB[1]-gV[1]*gB[0];
246
247 P[0]=gP[0]+gV[0]*sl+Ac*DVx;
248 P[1]=gP[1]+gV[1]*sl+Ac*DVy;
249 P[2]=gP[2]+gV[2]*sl+Ac*DVz;
250 V[0]=gV[0]+Av*DVx;
251 V[1]=gV[1]+Av*DVy;
252 V[2]=gV[2]+Av*DVz;
253
255
256 for(i=0;
i<3;
i++) gBf[i]=gB[i];
258 {
259 dBds[
i]=(gBf[
i]-gBi[
i])/sl;
261 }
262 nStep=nStepMax;
263 while(nStep>0)
264 {
265 c=D[0]*gP[0]+D[1]*gP[1]+D[2]*gP[2]+D[3];
266 b=D[0]*gV[0]+D[1]*gV[1]+D[2]*gV[2];
267 a=0.5*CQ*(gB[0]*(D[1]*gV[2]-D[2]*gV[1])+
268 gB[1]*(D[2]*gV[0]-D[0]*gV[2])+
269 gB[2]*(D[0]*gV[1]-D[1]*gV[0]));
270
272 if(fabs(ratio)>0.1)
273 useExpansion = false;
274 else useExpansion = true;
275
276 if(useExpansion) {
279 }
280 else {
282 if(descr<0.0)
283 {
284
285 return nullptr;
286 }
287 int signb = (
b<0.0)?-1:1;
288 sl = (-
b+signb*sqrt(descr))/(2*
a);
289 }
290
294 DVx=gV[1]*gB[2]-gV[2]*gB[1];
295 DVy=gV[2]*gB[0]-gV[0]*gB[2];
296 DVz=gV[0]*gB[1]-gV[1]*gB[0];
297
298 P[0]=gP[0]+gV[0]*
ds+Ac*DVx;
299 P[1]=gP[1]+gV[1]*
ds+Ac*DVy;
300 P[2]=gP[2]+gV[2]*
ds+Ac*DVz;
301 V[0]=gV[0]+Av*DVx;
302 V[1]=gV[1]+Av*DVy;
303 V[2]=gV[2]+Av*DVz;
305 {
306 gV[
i]=V[
i];gP[
i]=
P[
i];
307 }
308 for(i=0;
i<3;
i++) gB[i]+=dBds[i]*ds;
309 nStep--;
310 }
312 Rf[0]=lP[0];Rf[1]=lP[1];
313 Rf[2]=atan2(V[1],V[0]);
314
315 if(fabs(V[2])>1.0)
316 {
317 return nullptr;
318 }
319
320 Rf[3]=acos(V[2]);
322
323 gV[0]=sint*cosf;gV[1]=sint*sinf;gV[2]=
cost;
324
325 for(i=0;
i<4;
i++) D[i]=pSE->
getPar(i);
326 for(i=0;
i<3;
i++) gP[i]=gPi[i];
327
329 {
330 gB[
i]=0.5*(gBi[
i]+gBf[
i]);
331 }
332
333 c=D[0]*gP[0]+D[1]*gP[1]+D[2]*gP[2]+D[3];
334 b=D[0]*gV[0]+D[1]*gV[1]+D[2]*gV[2];
335 a=CQ*(gB[0]*(D[1]*gV[2]-D[2]*gV[1])+
336 gB[1]*(D[2]*gV[0]-D[0]*gV[2])+
337 gB[2]*(D[0]*gV[1]-D[1]*gV[0]));
338
340 if(fabs(ratio)>0.1)
341 useExpansion = false;
342 else useExpansion = true;
343
344 if(useExpansion) {
347 }
348 else {
350 if(descr<0.0)
351 {
352
353 return nullptr;
354 }
355 int signb = (
b<0.0)?-1:1;
356 s = (-
b+signb*sqrt(descr))/(2*
a);
357 }
358
362
363 DVx=gV[1]*gB[2]-gV[2]*gB[1];
364 DVy=gV[2]*gB[0]-gV[0]*gB[2];
365 DVz=gV[0]*gB[1]-gV[1]*gB[0];
366
367 P[0]=gP[0]+gV[0]*
s+Ac*DVx;
368 P[1]=gP[1]+gV[1]*
s+Ac*DVy;
369 P[2]=gP[2]+gV[2]*
s+Ac*DVz;
370
371 V[0]=gV[0]+Av*DVx;V[1]=gV[1]+Av*DVy;V[2]=gV[2]+Av*DVz;
372 if (std::abs(V[2]) > 1.0) {
373 return nullptr;
374 }
375
377
378 memset(&Jm[0][0],0,sizeof(Jm));
379
381
382 double coeff[3], dadVx,dadVy,dadVz,dadQ,dsdx,dsdy,dsdz,dsdVx,dsdVy,dsdVz,dsdQ;
383 coeff[0]=-
c*
c/(
b*
b*
b);
384 coeff[1]=
c*(1.0+3.0*
c*
a/(
b*
b))/(b*b);
385 coeff[2]=-(1.0+2.0*
c*
a/(
b*
b))/b;
386
387 dadVx=0.5*CQ*(-D[1]*gB[2]+D[2]*gB[1]);
388 dadVy=0.5*CQ*( D[0]*gB[2]-D[2]*gB[0]);
389 dadVz=0.5*CQ*(-D[0]*gB[1]+D[1]*gB[0]);
390 dadQ=0.5*
C*(D[0]*DVx+D[1]*DVy+D[2]*DVz);
391
392 dsdx=coeff[2]*D[0];
393 dsdy=coeff[2]*D[1];
394 dsdz=coeff[2]*D[2];
395 dsdVx=coeff[0]*dadVx+coeff[1]*D[0];
396 dsdVy=coeff[0]*dadVy+coeff[1]*D[1];
397 dsdVz=coeff[0]*dadVz+coeff[1]*D[2];
398 dsdQ=coeff[0]*dadQ;
399
400 Jm[0][0]=1.0+V[0]*dsdx;
401 Jm[0][1]= V[0]*dsdy;
402 Jm[0][2]= V[0]*dsdz;
403
404 Jm[0][3]=
s+V[0]*dsdVx;
405 Jm[0][4]= V[0]*dsdVy+Ac*gB[2];
406 Jm[0][5]= V[0]*dsdVz-Ac*gB[1];
407 Jm[0][6]= V[0]*dsdQ+Cc*DVx;
408
409 Jm[1][0]= V[1]*dsdx;
410 Jm[1][1]=1.0+V[1]*dsdy;
411 Jm[1][2]= V[1]*dsdz;
412
413 Jm[1][3]= V[1]*dsdVx-Ac*gB[2];
414 Jm[1][4]=
s+V[1]*dsdVy;
415 Jm[1][5]= V[1]*dsdVz+Ac*gB[0];
416 Jm[1][6]= V[1]*dsdQ+Cc*DVy;
417
418 Jm[2][0]= V[2]*dsdx;
419 Jm[2][1]= V[2]*dsdy;
420 Jm[2][2]=1.0+V[2]*dsdz;
421 Jm[2][3]= V[2]*dsdVx+Ac*gB[1];
422 Jm[2][4]= V[2]*dsdVy-Ac*gB[0];
423 Jm[2][5]=
s+V[2]*dsdVz;
424 Jm[2][6]= V[2]*dsdQ+Cc*DVz;
425
426 Jm[3][0]=dsdx*CQ*DVx;
427 Jm[3][1]=dsdy*CQ*DVx;
428 Jm[3][2]=dsdz*CQ*DVx;
429
430 Jm[3][3]=1.0+dsdVx*CQ*DVx;
431 Jm[3][4]=CQ*(dsdVy*DVx+
s*gB[2]);
432 Jm[3][5]=CQ*(dsdVz*DVx-
s*gB[1]);
433
434 Jm[3][6]=(CQ*dsdQ+
C*
s)*DVx;
435
436 Jm[4][0]=dsdx*CQ*DVy;
437 Jm[4][1]=dsdy*CQ*DVy;
438 Jm[4][2]=dsdz*CQ*DVy;
439
440 Jm[4][3]=CQ*(dsdVx*DVy-
s*gB[2]);
441 Jm[4][4]=1.0+dsdVy*CQ*DVy;
442 Jm[4][5]=CQ*(dsdVz*DVy+
s*gB[0]);
443
444 Jm[4][6]=(CQ*dsdQ+
C*
s)*DVy;
445
446 Jm[5][0]=dsdx*CQ*DVz;
447 Jm[5][1]=dsdy*CQ*DVz;
448 Jm[5][2]=dsdz*CQ*DVz;
449 Jm[5][3]=CQ*(dsdVx*DVz+
s*gB[1]);
450 Jm[5][4]=CQ*(dsdVy*DVz-
s*gB[0]);
451 Jm[5][5]=1.0+dsdVz*CQ*DVz;
452 Jm[5][6]=(CQ*dsdQ+
C*
s)*DVz;
453
454 Jm[6][6]=1.0;
455
456 memset(&J1[0][0],0,sizeof(J1));
457
458 J1[0][0]=M[0][0];J1[0][1]=M[0][1];J1[0][2]=M[0][2];
459 J1[1][0]=M[1][0];J1[1][1]=M[1][1];J1[1][2]=M[1][2];
460 J1[2][3]=-V[1]/(V[0]*V[0]+V[1]*V[1]);
461 J1[2][4]= V[0]/(V[0]*V[0]+V[1]*V[1]);
462 J1[3][5]=-1.0/sqrt(1-V[2]*V[2]);
463 J1[4][6]=1.0;
464
466 {
468 Buf[j][i]=J1[j][0]*Jm[0][i]+J1[j][1]*Jm[1][i]+J1[j][2]*Jm[2][i];
469 Buf[2][
i]=J1[2][3]*Jm[3][
i]+J1[2][4]*Jm[4][
i];
470 Buf[3][
i]=J1[3][5]*Jm[5][
i];
472 }
473
474 if(pSB!=nullptr)
475 {
477 {
478 J[
i][0]=Buf[
i][0]*J0[0][0]+Buf[
i][1]*J0[1][0]+Buf[
i][2]*J0[2][0];
479 J[
i][1]=Buf[
i][0]*J0[0][1]+Buf[
i][1]*J0[1][1]+Buf[
i][2]*J0[2][1];
480 J[
i][2]=Buf[
i][3]*J0[3][2]+Buf[
i][4]*J0[4][2];
481 J[
i][3]=Buf[
i][3]*J0[3][3]+Buf[
i][4]*J0[4][3]+Buf[
i][5]*J0[5][3];
483 }
484 }
485 else
486 {
488 {
489 J[
i][0]=Buf[
i][0]*J0[0][0]+Buf[
i][1]*J0[1][0];
491 J[
i][2]=Buf[
i][0]*J0[0][2]+Buf[
i][1]*J0[1][2]+Buf[
i][3]*J0[3][2]+Buf[
i][4]*J0[4][2];
492 J[
i][3]=Buf[
i][3]*J0[3][3]+Buf[
i][4]*J0[4][3]+Buf[
i][5]*J0[5][3];
494 }
495 }
496 }
497 else {
503 memset(&J[0][0],0,sizeof(J));
504 for(i=0;
i<5;
i++) J[i][i]=1.0;
505 }
506
507 for(i=0;
i<5;
i++)
for(j=0;
j<5;
j++)
508 {
510 }
511 for(i=0;
i<5;
i++)
for(j=i;
j<5;
j++)
512 {
514 for(m=0;
m<5;
m++) Gf[i][j]+=AG[i][m]*J[j][m];
516 }
517
518 Trk::TrkTrackState* pTE=new Trk::TrkTrackState(pTS);
519
520
521 double Rtmp[5];
522
523 for(i=0;
i<4;
i++) Rtmp[i] = Rf[i];
524 Rtmp[4] = 0.001*Rf[4];
525
529
532
534
537
538
539
540 for(
int idx=0;
idx<5;
idx++) {
543 delete pTE;
544 return nullptr;
545 }
546 }
547
549 for(i=0;i<5;i++) for(j=i;j<5;j++)
550 {
552 }
553 Gi = Gi.inverse();
554
555 for(i=0;
i<5;
i++)
for(j=0;
j<5;
j++)
556 {
558 for(m=0;
m<5;
m++) A[i][j]+=AG[m][i]*Gi(m,j);
559 }
562
563 return pTE;
564}
#define ATH_MSG_DEBUG(x,...)
void diff(const Jet &rJet1, const Jet &rJet2, std::map< std::string, double > varDiff)
Difference between jets - Non-Class function required by trigger.
void getMagneticField(double[3], double *, MagField::AtlasFieldCache &) const
void transformPointToLocal(const double *, double *)
double getRotMatrix(int, int)
double getInvRotMatrix(int, int)
void transformPointToGlobal(const double *, double *)
void applyEnergyLoss(int)
void setPreviousState(TrkTrackState *)
void applyMultipleScattering()
void setSmootherGain(double A[5][5])
void attachToSurface(TrkPlanarSurface *)
int cost(std::vector< std::string > &files, node &n, const std::string &directory="", bool deleteref=false, bool relocate=false)