137{
138 const double C=0.02999975/1000.0;
139 const double minStep=30.0;
140
141 double J[5][5],Rf[5],AG[5][5],Gf[5][5],
A[5][5];
143
144 bool samePlane=false;
145
146 if(pSB!=nullptr)
147 {
151 samePlane=true;
153 }
154 }
155
156 if(!samePlane) {
157
158 double gP[3],gPi[3],lP[3],gV[3],
a,
b,
c,
s,J0[7][5],
descr,CQ,Ac,Av,Cc;
159 double V[3],
P[3],M[3][3],D[4],Jm[7][7],
160 J1[5][7],gB[3],gBi[3],gBf[3],dBds[3],Buf[5][7],DVx,DVy,DVz;
161 int nStep,nStepMax;
163
164 double sint,
cost,sinf,cosf;
168
169 memset(&J0[0][0],0,sizeof(J0));
170
171 if(pSB!=nullptr)
172 {
177
178 J0[0][0]=
L[0][0];J0[0][1]=
L[0][1];
179 J0[1][0]=
L[1][0];J0[1][1]=
L[1][1];
180 J0[2][0]=
L[2][0];J0[2][1]=
L[2][1];
181 J0[3][2]=-sinf*sint;J0[3][3]=cosf*
cost;
182 J0[4][2]= cosf*sint;J0[4][3]=sinf*
cost;
183 J0[5][3]=-sint;
184 J0[6][4]=1.0;
185 }
186 else
187 {
193 J0[2][1]=1.0;
194 J0[3][2]=-sinf*sint;J0[3][3]=cosf*
cost;
195 J0[4][2]= cosf*sint;J0[4][3]=sinf*
cost;
196 J0[5][3]=-sint;
197 J0[6][4]=1.0;
198 }
199 for(i=0;
i<4;
i++) D[i]=pSE->
getPar(i);
200 for(i=0;
i<3;
i++) gPi[i]=gP[i];
201
203
204 for(i=0;
i<3;
i++) gBi[i]=gB[i];
205
206 c=D[0]*gP[0]+D[1]*gP[1]+D[2]*gP[2]+D[3];
207 b=D[0]*gV[0]+D[1]*gV[1]+D[2]*gV[2];
208 a=0.5*CQ*(gB[0]*(D[1]*gV[2]-D[2]*gV[1])+
209 gB[1]*(D[2]*gV[0]-D[0]*gV[2])+
210 gB[2]*(D[0]*gV[1]-D[1]*gV[0]));
211
213
214 if(descr<0.0)
215 {
216
217 return nullptr;
218 }
219
220 bool useExpansion=true;
222
223 if(fabs(ratio)>0.1)
224 useExpansion = false;
225
226 if(useExpansion) {
229 }
230 else {
231 int signb = (
b<0.0)?-1:1;
232 sl = (-
b+signb*sqrt(descr))/(2*
a);
233 }
234
235 if(fabs(sl)<minStep) nStepMax=1;
236 else
237 {
238 nStepMax=(
int)(fabs(sl)/minStep)+1;
239 }
240 if((nStepMax<0)||(nStepMax>1000))
241 {
242 return nullptr;
243 }
244 Av=sl*CQ;
245 Ac=0.5*sl*Av;
246 DVx=gV[1]*gB[2]-gV[2]*gB[1];
247 DVy=gV[2]*gB[0]-gV[0]*gB[2];
248 DVz=gV[0]*gB[1]-gV[1]*gB[0];
249
250 P[0]=gP[0]+gV[0]*sl+Ac*DVx;
251 P[1]=gP[1]+gV[1]*sl+Ac*DVy;
252 P[2]=gP[2]+gV[2]*sl+Ac*DVz;
253 V[0]=gV[0]+Av*DVx;
254 V[1]=gV[1]+Av*DVy;
255 V[2]=gV[2]+Av*DVz;
256
258
259 for(i=0;
i<3;
i++) gBf[i]=gB[i];
261 {
262 dBds[
i]=(gBf[
i]-gBi[
i])/sl;
264 }
265 nStep=nStepMax;
266 while(nStep>0)
267 {
268 c=D[0]*gP[0]+D[1]*gP[1]+D[2]*gP[2]+D[3];
269 b=D[0]*gV[0]+D[1]*gV[1]+D[2]*gV[2];
270 a=0.5*CQ*(gB[0]*(D[1]*gV[2]-D[2]*gV[1])+
271 gB[1]*(D[2]*gV[0]-D[0]*gV[2])+
272 gB[2]*(D[0]*gV[1]-D[1]*gV[0]));
273
275 if(fabs(ratio)>0.1)
276 useExpansion = false;
277 else useExpansion = true;
278
279 if(useExpansion) {
282 }
283 else {
285 if(descr<0.0)
286 {
287
288 return nullptr;
289 }
290 int signb = (
b<0.0)?-1:1;
291 sl = (-
b+signb*sqrt(descr))/(2*
a);
292 }
293
297 DVx=gV[1]*gB[2]-gV[2]*gB[1];
298 DVy=gV[2]*gB[0]-gV[0]*gB[2];
299 DVz=gV[0]*gB[1]-gV[1]*gB[0];
300
301 P[0]=gP[0]+gV[0]*
ds+Ac*DVx;
302 P[1]=gP[1]+gV[1]*
ds+Ac*DVy;
303 P[2]=gP[2]+gV[2]*
ds+Ac*DVz;
304 V[0]=gV[0]+Av*DVx;
305 V[1]=gV[1]+Av*DVy;
306 V[2]=gV[2]+Av*DVz;
308 {
309 gV[
i]=V[
i];gP[
i]=
P[
i];
310 }
311 for(i=0;
i<3;
i++) gB[i]+=dBds[i]*ds;
312 nStep--;
313 }
315 Rf[0]=lP[0];Rf[1]=lP[1];
316 Rf[2]=atan2(V[1],V[0]);
317
318 if(fabs(V[2])>1.0)
319 {
320 return nullptr;
321 }
322
323 Rf[3]=acos(V[2]);
325
326 gV[0]=sint*cosf;gV[1]=sint*sinf;gV[2]=
cost;
327
328 for(i=0;
i<4;
i++) D[i]=pSE->
getPar(i);
329 for(i=0;
i<3;
i++) gP[i]=gPi[i];
330
332 {
333 gB[
i]=0.5*(gBi[
i]+gBf[
i]);
334 }
335
336 c=D[0]*gP[0]+D[1]*gP[1]+D[2]*gP[2]+D[3];
337 b=D[0]*gV[0]+D[1]*gV[1]+D[2]*gV[2];
338 a=CQ*(gB[0]*(D[1]*gV[2]-D[2]*gV[1])+
339 gB[1]*(D[2]*gV[0]-D[0]*gV[2])+
340 gB[2]*(D[0]*gV[1]-D[1]*gV[0]));
341
343 if(fabs(ratio)>0.1)
344 useExpansion = false;
345 else useExpansion = true;
346
347 if(useExpansion) {
350 }
351 else {
353 if(descr<0.0)
354 {
355
356 return nullptr;
357 }
358 int signb = (
b<0.0)?-1:1;
359 s = (-
b+signb*sqrt(descr))/(2*
a);
360 }
361
365
366 DVx=gV[1]*gB[2]-gV[2]*gB[1];
367 DVy=gV[2]*gB[0]-gV[0]*gB[2];
368 DVz=gV[0]*gB[1]-gV[1]*gB[0];
369
370 P[0]=gP[0]+gV[0]*
s+Ac*DVx;
371 P[1]=gP[1]+gV[1]*
s+Ac*DVy;
372 P[2]=gP[2]+gV[2]*
s+Ac*DVz;
373
374 V[0]=gV[0]+Av*DVx;V[1]=gV[1]+Av*DVy;V[2]=gV[2]+Av*DVz;
375 if (std::abs(V[2]) > 1.0) {
376 return nullptr;
377 }
378
380
381 memset(&Jm[0][0],0,sizeof(Jm));
382
384
385 double coeff[3], dadVx,dadVy,dadVz,dadQ,dsdx,dsdy,dsdz,dsdVx,dsdVy,dsdVz,dsdQ;
386 coeff[0]=-
c*
c/(
b*
b*
b);
387 coeff[1]=
c*(1.0+3.0*
c*
a/(
b*
b))/(b*b);
388 coeff[2]=-(1.0+2.0*
c*
a/(
b*
b))/b;
389
390 dadVx=0.5*CQ*(-D[1]*gB[2]+D[2]*gB[1]);
391 dadVy=0.5*CQ*( D[0]*gB[2]-D[2]*gB[0]);
392 dadVz=0.5*CQ*(-D[0]*gB[1]+D[1]*gB[0]);
393 dadQ=0.5*
C*(D[0]*DVx+D[1]*DVy+D[2]*DVz);
394
395 dsdx=coeff[2]*D[0];
396 dsdy=coeff[2]*D[1];
397 dsdz=coeff[2]*D[2];
398 dsdVx=coeff[0]*dadVx+coeff[1]*D[0];
399 dsdVy=coeff[0]*dadVy+coeff[1]*D[1];
400 dsdVz=coeff[0]*dadVz+coeff[1]*D[2];
401 dsdQ=coeff[0]*dadQ;
402
403 Jm[0][0]=1.0+V[0]*dsdx;
404 Jm[0][1]= V[0]*dsdy;
405 Jm[0][2]= V[0]*dsdz;
406
407 Jm[0][3]=
s+V[0]*dsdVx;
408 Jm[0][4]= V[0]*dsdVy+Ac*gB[2];
409 Jm[0][5]= V[0]*dsdVz-Ac*gB[1];
410 Jm[0][6]= V[0]*dsdQ+Cc*DVx;
411
412 Jm[1][0]= V[1]*dsdx;
413 Jm[1][1]=1.0+V[1]*dsdy;
414 Jm[1][2]= V[1]*dsdz;
415
416 Jm[1][3]= V[1]*dsdVx-Ac*gB[2];
417 Jm[1][4]=
s+V[1]*dsdVy;
418 Jm[1][5]= V[1]*dsdVz+Ac*gB[0];
419 Jm[1][6]= V[1]*dsdQ+Cc*DVy;
420
421 Jm[2][0]= V[2]*dsdx;
422 Jm[2][1]= V[2]*dsdy;
423 Jm[2][2]=1.0+V[2]*dsdz;
424 Jm[2][3]= V[2]*dsdVx+Ac*gB[1];
425 Jm[2][4]= V[2]*dsdVy-Ac*gB[0];
426 Jm[2][5]=
s+V[2]*dsdVz;
427 Jm[2][6]= V[2]*dsdQ+Cc*DVz;
428
429 Jm[3][0]=dsdx*CQ*DVx;
430 Jm[3][1]=dsdy*CQ*DVx;
431 Jm[3][2]=dsdz*CQ*DVx;
432
433 Jm[3][3]=1.0+dsdVx*CQ*DVx;
434 Jm[3][4]=CQ*(dsdVy*DVx+
s*gB[2]);
435 Jm[3][5]=CQ*(dsdVz*DVx-
s*gB[1]);
436
437 Jm[3][6]=(CQ*dsdQ+
C*
s)*DVx;
438
439 Jm[4][0]=dsdx*CQ*DVy;
440 Jm[4][1]=dsdy*CQ*DVy;
441 Jm[4][2]=dsdz*CQ*DVy;
442
443 Jm[4][3]=CQ*(dsdVx*DVy-
s*gB[2]);
444 Jm[4][4]=1.0+dsdVy*CQ*DVy;
445 Jm[4][5]=CQ*(dsdVz*DVy+
s*gB[0]);
446
447 Jm[4][6]=(CQ*dsdQ+
C*
s)*DVy;
448
449 Jm[5][0]=dsdx*CQ*DVz;
450 Jm[5][1]=dsdy*CQ*DVz;
451 Jm[5][2]=dsdz*CQ*DVz;
452 Jm[5][3]=CQ*(dsdVx*DVz+
s*gB[1]);
453 Jm[5][4]=CQ*(dsdVy*DVz-
s*gB[0]);
454 Jm[5][5]=1.0+dsdVz*CQ*DVz;
455 Jm[5][6]=(CQ*dsdQ+
C*
s)*DVz;
456
457 Jm[6][6]=1.0;
458
459 memset(&J1[0][0],0,sizeof(J1));
460
461 J1[0][0]=M[0][0];J1[0][1]=M[0][1];J1[0][2]=M[0][2];
462 J1[1][0]=M[1][0];J1[1][1]=M[1][1];J1[1][2]=M[1][2];
463 J1[2][3]=-V[1]/(V[0]*V[0]+V[1]*V[1]);
464 J1[2][4]= V[0]/(V[0]*V[0]+V[1]*V[1]);
465 J1[3][5]=-1.0/sqrt(1-V[2]*V[2]);
466 J1[4][6]=1.0;
467
469 {
470 for(j=0;j<2;j++)
471 Buf[j][i]=J1[j][0]*Jm[0][i]+J1[j][1]*Jm[1][i]+J1[j][2]*Jm[2][i];
472 Buf[2][
i]=J1[2][3]*Jm[3][
i]+J1[2][4]*Jm[4][
i];
473 Buf[3][
i]=J1[3][5]*Jm[5][
i];
475 }
476
477 if(pSB!=nullptr)
478 {
480 {
481 J[
i][0]=Buf[
i][0]*J0[0][0]+Buf[
i][1]*J0[1][0]+Buf[
i][2]*J0[2][0];
482 J[
i][1]=Buf[
i][0]*J0[0][1]+Buf[
i][1]*J0[1][1]+Buf[
i][2]*J0[2][1];
483 J[
i][2]=Buf[
i][3]*J0[3][2]+Buf[
i][4]*J0[4][2];
484 J[
i][3]=Buf[
i][3]*J0[3][3]+Buf[
i][4]*J0[4][3]+Buf[
i][5]*J0[5][3];
486 }
487 }
488 else
489 {
491 {
492 J[
i][0]=Buf[
i][0]*J0[0][0]+Buf[
i][1]*J0[1][0];
494 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];
495 J[
i][3]=Buf[
i][3]*J0[3][3]+Buf[
i][4]*J0[4][3]+Buf[
i][5]*J0[5][3];
497 }
498 }
499 }
500 else {
506 memset(&J[0][0],0,sizeof(J));
507 for(i=0;
i<5;
i++) J[i][i]=1.0;
508 }
509
510 for(i=0;
i<5;
i++)
for(j=0;j<5;j++)
511 {
513 }
514 for(i=0;
i<5;
i++)
for(j=i;j<5;j++)
515 {
517 for(m=0;
m<5;
m++) Gf[i][j]+=AG[i][m]*J[j][m];
519 }
520
521 Trk::TrkTrackState* pTE=new Trk::TrkTrackState(pTS);
522
523
524 double Rtmp[5];
525
526 for(i=0;
i<4;
i++) Rtmp[i] = Rf[i];
527 Rtmp[4] = 0.001*Rf[4];
528
532
535
537
540
541
542
543 for(
int idx=0;
idx<5;
idx++) {
546 delete pTE;
547 return nullptr;
548 }
549 }
550
552 for(i=0;i<5;i++) for(j=i;j<5;j++)
553 {
555 }
556 Gi = Gi.inverse();
557
558 for(i=0;
i<5;
i++)
for(j=0;j<5;j++)
559 {
561 for(m=0;
m<5;
m++) A[i][j]+=AG[m][i]*Gi(m,j);
562 }
565
566 return pTE;
567}
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 *)
void transformPointToGlobal(const double *, double *)
double getInvRotMatrix(int, int)
double getRotMatrix(int, int)
void applyMultipleScattering()
void setPreviousState(TrkTrackState *)
void attachToSurface(TrkPlanarSurface *)
void applyEnergyLoss(int)
void setSmootherGain(double A[5][5])
int cost(std::vector< std::string > &files, node &n, const std::string &directory="", bool deleteref=false, bool relocate=false)