70 {
71
72 double xc, yc;
73 double rho;
74 double xx[2], yy[2];
75 double tofx, tofy;
76 double nx, ny;
77 double cosx, sinx;
78 double aa, ca, cb, cc, detm;
79
80
81 int Id_cdfztrk = id;
83
84 if ( Id_cdfztrk < 0 )
85 {
86 if ( _debug )
87 std::cout << " TofFz_helix *** wrong id_cdfztrk =" << Id_cdfztrk << std::endl;
88 return ( -1 );
89 }
90
92 HepVector a( 5, 0 );
93 a[0] = fzisan[0];
94 a[1] = fzisan[1];
95 a[2] = fzisan[2];
96 a[3] = fzisan[3];
97 a[4] = fzisan[4];
98
99
100 if (
abs( a[3] ) > 50.0 ||
abs( a[4] ) > 500.0 )
return ( -5 );
101
102
103
104 Helix helix1( pivot1, a );
106 Phi0 = helix1.a()[1];
107 Kappa = helix1.a()[2];
109 Tanl = helix1.a()[4];
110
111
112 rho = helix1.radius();
113
115 xc = xyz( 0 );
116 yc = xyz( 1 );
117
118
119 Hep3Vector vxyz = helix1.direction();
120 nx = vxyz( 0 );
121 ny = vxyz( 1 );
122
123
124 if ( fabs( rho ) >
R_tof / 2 )
125 {
126
127
128
129 if ( xc == 0.0 && yc == 0.0 ) return ( -3 );
130
131 ca = xc * xc + yc * yc;
132 aa = 0.5 * ( ca - rho * rho +
R_tof *
R_tof );
133
134 if ( xc != 0.0 )
135 {
136 cb = aa * yc;
138
139 detm = cb * cb - ca * cc;
140 if ( detm > 0.0 )
141 {
142 yy[0] = ( cb + sqrt( detm ) ) / ca;
143 xx[0] = ( aa - yy[0] * yc ) / xc;
144 yy[1] = ( cb - sqrt( detm ) ) / ca;
145 xx[1] = ( aa - yy[1] * yc ) / xc;
146 }
147 else return ( -1 );
148 }
149 else
150 {
151 cb = aa * xc;
153
154 detm = cb * cb - ca * cc;
155 if ( detm > 0.0 )
156 {
157 xx[0] = ( cb + sqrt( detm ) ) / ca;
158 yy[0] = ( aa - xx[0] * xc ) / yc;
159 xx[1] = ( cb - sqrt( detm ) ) / ca;
160 yy[1] = ( aa - xx[1] * xc ) / yc;
161 }
162 else return ( -2 );
163 }
164
165
166
167 if ( xx[0] * nx + yy[0] * ny > 0.0 )
168 {
169 tofx = xx[0];
170 tofy = yy[0];
171 }
172 else
173 {
174 tofx = xx[1];
175 tofy = yy[1];
176 }
177
178 double fi = atan2( tofy, tofx );
179
180 if ( fi < 0.0 ) fi = pi2 + fi;
182
184
185
187 helix1.pivot( pivot2 );
188
189
190 Phi1 = helix1.a()[1];
191 Z_tof = helix1.a()[3];
192
193 double dphi =
194 ( -xc * ( tofx - xc ) - yc * ( tofy - yc ) ) / sqrt( xc * xc + yc * yc ) / fabs( rho );
195 dphi = acos( dphi );
196 Pathl = fabs( rho * dphi );
197 double coscor = sqrt( 1.0 +
Tanl *
Tanl );
199
200
201 Hep3Vector vxyz1 = helix1.direction();
202 nx = vxyz1( 0 );
203 ny = vxyz1( 1 );
204
205 double corxy = ( nx * tofx + ny * tofy ) /
R_tof;
207
208 if (
abs(
Z_tof ) < 115 )
return ( 1 );
209
210
211
212
214 {
215
222
223 double x_endtof =
225 double y_endtof =
227 r_endtof = sqrt( x_endtof * x_endtof + y_endtof * y_endtof );
228
229 Helix helix1( pivot1, a );
230
231 double x1_endtof = helix1.x( phi ).x();
232 double y1_endtof = helix1.x( phi ).y();
233 double z1_endtof = helix1.x( phi ).z();
234
235 double fi_endtof = atan2( y_endtof, x_endtof );
236 if ( fi_endtof < 0 ) fi_endtof = pi2 + fi_endtof;
237
238 Tofid = (int)( fi_endtof / piby24 );
240
241
242
243 double Z_etf_1 = 133.673;
244 double Pathl_1 = Z_etf_1 /
sin( atan(
Tanl ) );
245 double phi0_1 = -( Z_etf_1 /
Tanl ) / rho;
246 double phi_1 = -( Z_etf_1 -
Dz ) /
Tanl / rho;
247 double phi1_1 = ( Z_etf_1 *
Kappa * 0.003 ) /
Tanl;
248 double phi2_1 = ( Z_etf_1 -
Dz ) *
Kappa * 0.003 /
Tanl;
249
250 double x1_etf =
252 double y1_etf =
254 double r_etf_1 = sqrt( x1_etf * x1_etf + y1_etf * y1_etf );
255
256 double x_etf_1 = helix1.x( phi_1 ).x();
257 double y_etf_1 = helix1.x( phi_1 ).y();
258 double z_etf_1 = helix1.x( phi_1 ).z();
259
260
261 double fi_etf_1 = atan2( y1_etf, x1_etf ) + 1.25 * piby1 / 180.0;
262 if ( fi_etf_1 < 0 ) { fi_etf_1 = pi2 + fi_etf_1; }
263
264 int Etfid_1 = (int)( fi_etf_1 / piby18 );
265 if ( Etfid_1 > 35 ) { Etfid_1 = Etfid_1 - 36; }
266
267
268 if ( Etfid_1 % 2 == 1 )
269 {
274 }
275
276 else
277 {
278
279 double Z_etf_2 = 136.573;
280 double Pathl_2 = Z_etf_2 /
sin( atan(
Tanl ) );
281 double phi0_2 = -( Z_etf_2 /
Tanl ) / rho;
282 double phi_2 = -( Z_etf_2 -
Dz ) /
Tanl / rho;
283 double phi1_2 = ( Z_etf_2 *
Kappa * 0.003 ) /
Tanl;
284 double phi2_2 = ( Z_etf_2 -
Dz ) *
Kappa * 0.003 /
Tanl;
285
290 double r_etf_2 = sqrt( x2_etf * x2_etf + y2_etf * y2_etf );
291
292 double x_etf_2 = helix1.x( phi_2 ).x();
293 double y_etf_2 = helix1.x( phi_2 ).y();
294 double z_etf_2 = helix1.x( phi_2 ).z();
295
296
297 double fi_etf_2 = atan2( y2_etf, x2_etf ) + 1.25 * piby1 / 180.0;
298 if ( fi_etf_2 < 0 ) { fi_etf_2 = pi2 + fi_etf_2; }
299
300 int Etfid_2 = (int)( fi_etf_2 / piby18 );
301 if ( Etfid_2 > 35 ) { Etfid_2 = Etfid_2 - 36; }
302
303 if ( Etfid_2 % 2 == 1 )
304 {
305 int tmp1 = (int)( ( fi_etf_2 + 0.5 * piby18 ) / piby18 );
306 int tmp2 = (int)( ( fi_etf_2 - 0.5 * piby18 ) / piby18 );
307 if ( tmp1 == Etfid_2 ) { Etfid_2 = tmp2; }
308 else { Etfid_2 = tmp1; }
309 }
310
315 }
316
317 return ( 0 );
318 }
319 else if (
Z_tof <= -115 )
320 {
321
328
329 double x_endtof =
331 double y_endtof =
333 r_endtof = sqrt( x_endtof * x_endtof + y_endtof * y_endtof );
334
335 Helix helix1( pivot1, a );
336
337 double x1_endtof = helix1.x( phi ).x();
338 double y1_endtof = helix1.x( phi ).y();
339 double z1_endtof = helix1.x( phi ).z();
340
341 double fi_endtof = atan2( y_endtof, x_endtof );
342 if ( fi_endtof < 0 ) fi_endtof = pi2 + fi_endtof;
343
344 Tofid = (int)( fi_endtof / piby24 );
346
347
348 double Z_etf_1 = -133.673;
349 double Pathl_1 = Z_etf_1 /
sin( atan(
Tanl ) );
350 double phi0_1 = -( Z_etf_1 /
Tanl ) / rho;
351 double phi_1 = -( Z_etf_1 -
Dz ) /
Tanl / rho;
352 double phi1_1 = ( Z_etf_1 *
Kappa * 0.003 ) /
Tanl;
353 double phi2_1 = ( Z_etf_1 -
Dz ) *
Kappa * 0.003 /
Tanl;
354
355 double x1_etf =
357 double y1_etf =
359 double r_etf_1 = sqrt( x1_etf * x1_etf + y1_etf * y1_etf );
360
361 double x_etf_1 = helix1.x( phi_1 ).x();
362 double y_etf_1 = helix1.x( phi_1 ).y();
363 double z_etf_1 = helix1.x( phi_1 ).z();
364
365
366 double fi_etf_1 = atan2( y1_etf, x1_etf ) - 11.25 * piby1 / 180.0;
367 if ( fi_etf_1 < 0 ) { fi_etf_1 = pi2 + fi_etf_1; }
368
369 int Etfid_1 = (int)( fi_etf_1 / piby18 );
370 if ( Etfid_1 > 35 ) { Etfid_1 = Etfid_1 - 36; }
371
372
373 if ( Etfid_1 % 2 == 1 )
374 {
379 }
380
381 else
382 {
383 double Z_etf_2 = -136.573;
384 double Pathl_2 = Z_etf_2 /
sin( atan(
Tanl ) );
385 double phi0_2 = -( Z_etf_2 /
Tanl ) / rho;
386 double phi_2 = -( Z_etf_2 -
Dz ) /
Tanl / rho;
387 double phi1_2 = ( Z_etf_2 *
Kappa * 0.003 ) /
Tanl;
388 double phi2_2 = ( Z_etf_2 -
Dz ) *
Kappa * 0.003 /
Tanl;
389
394 double r_etf_2 = sqrt( x2_etf * x2_etf + y2_etf * y2_etf );
395
396 double x_etf_2 = helix1.x( phi_2 ).x();
397 double y_etf_2 = helix1.x( phi_2 ).y();
398 double z_etf_2 = helix1.x( phi_2 ).z();
399
400
401 double fi_etf_2 = atan2( y2_etf, x2_etf ) - 11.25 * piby1 / 180.0;
402 if ( fi_etf_2 < 0 ) { fi_etf_2 = pi2 + fi_etf_2; }
403
404 int Etfid_2 = (int)( fi_etf_2 / piby18 );
405 if ( Etfid_2 > 35 ) { Etfid_2 = Etfid_2 - 36; }
406
407 if ( Etfid_2 % 2 == 1 )
408 {
409 int tmp1 = (int)( ( fi_etf_2 + 0.5 * piby18 ) / piby18 );
410 int tmp2 = (int)( ( fi_etf_2 - 0.5 * piby18 ) / piby18 );
411 if ( tmp1 == Etfid_2 ) { Etfid_2 = tmp2; }
412 else { Etfid_2 = tmp1; }
413 }
414
419 }
420
421 return ( 2 );
422 }
423 }
424 else
425 {
427 {
428
435
436 double x_endtof =
438 double y_endtof =
440 r_endtof = sqrt( x_endtof * x_endtof + y_endtof * y_endtof );
441
442 double fi_endtof = atan2( y_endtof, x_endtof );
443 if ( fi_endtof < 0 ) fi_endtof = pi2 + fi_endtof;
444 Tofid = (int)( fi_endtof / piby24 );
446
447
448 double Z_etf_1 = 133.673;
449 double Pathl_1 = Z_etf_1 /
sin( atan(
Tanl ) );
450 double phi0_1 = -( Z_etf_1 /
Tanl ) / rho;
451 double phi_1 = -( Z_etf_1 -
Dz ) /
Tanl / rho;
452 double phi1_1 = ( Z_etf_1 *
Kappa * 0.003 ) /
Tanl;
453 double phi2_1 = ( Z_etf_1 -
Dz ) *
Kappa * 0.003 /
Tanl;
454
455 double x1_etf =
457 double y1_etf =
459 double r_etf_1 = sqrt( x1_etf * x1_etf + y1_etf * y1_etf );
460
461 double x_etf_1 = helix1.x( phi_1 ).x();
462 double y_etf_1 = helix1.x( phi_1 ).y();
463 double z_etf_1 = helix1.x( phi_1 ).z();
464
465
466 double fi_etf_1 = atan2( y1_etf, x1_etf ) + 1.25 * piby1 / 180.0;
467 if ( fi_etf_1 < 0 ) { fi_etf_1 = pi2 + fi_etf_1; }
468
469 int Etfid_1 = (int)( fi_etf_1 / piby18 );
470 if ( Etfid_1 > 35 ) { Etfid_1 = Etfid_1 - 36; }
471
472
473 if ( Etfid_1 % 2 == 1 )
474 {
479 }
480
481 else
482 {
483 double Z_etf_2 = 136.573;
484 double Pathl_2 = Z_etf_2 /
sin( atan(
Tanl ) );
485 double phi0_2 = -( Z_etf_2 /
Tanl ) / rho;
486 double phi_2 = -( Z_etf_2 -
Dz ) /
Tanl / rho;
487 double phi1_2 = ( Z_etf_2 *
Kappa * 0.003 ) /
Tanl;
488 double phi2_2 = ( Z_etf_2 -
Dz ) *
Kappa * 0.003 /
Tanl;
489
494 double r_etf_2 = sqrt( x2_etf * x2_etf + y2_etf * y2_etf );
495
496 double x_etf_2 = helix1.x( phi_2 ).x();
497 double y_etf_2 = helix1.x( phi_2 ).y();
498 double z_etf_2 = helix1.x( phi_2 ).z();
499
500
501 double fi_etf_2 = atan2( y2_etf, x2_etf ) + 1.25 * piby1 / 180.0;
502 if ( fi_etf_2 < 0 ) { fi_etf_2 = pi2 + fi_etf_2; }
503
504 int Etfid_2 = (int)( fi_etf_2 / piby18 );
505 if ( Etfid_2 > 35 ) { Etfid_2 = Etfid_2 - 36; }
506
507 if ( Etfid_2 % 2 == 1 )
508 {
509 int tmp1 = (int)( ( fi_etf_2 + 0.5 * piby18 ) / piby18 );
510 int tmp2 = (int)( ( fi_etf_2 - 0.5 * piby18 ) / piby18 );
511 if ( tmp1 == Etfid_2 ) { Etfid_2 = tmp2; }
512 else { Etfid_2 = tmp1; }
513 }
518 }
519
520 return ( 0 );
521 }
522 else
523 {
524
531
532 double x_endtof =
534 double y_endtof =
536 r_endtof = sqrt( x_endtof * x_endtof + y_endtof * y_endtof );
537
538 double fi_endtof = atan2( y_endtof, x_endtof );
539 if ( fi_endtof < 0 ) fi_endtof = pi2 + fi_endtof;
540 Tofid = (int)( fi_endtof / piby24 );
542
543
544 double Z_etf_1 = -133.673;
545 double Pathl_1 = Z_etf_1 /
sin( atan(
Tanl ) );
546 double phi0_1 = -( Z_etf_1 /
Tanl ) / rho;
547 double phi_1 = -( Z_etf_1 -
Dz ) /
Tanl / rho;
548 double phi1_1 = ( Z_etf_1 *
Kappa * 0.003 ) /
Tanl;
549 double phi2_1 = ( Z_etf_1 -
Dz ) *
Kappa * 0.003 /
Tanl;
550
551 double x1_etf =
553 double y1_etf =
555 double r_etf_1 = sqrt( x1_etf * x1_etf + y1_etf * y1_etf );
556
557 double x_etf_1 = helix1.x( phi_1 ).x();
558 double y_etf_1 = helix1.x( phi_1 ).y();
559 double z_etf_1 = helix1.x( phi_1 ).z();
560
561
562 double fi_etf_1 = atan2( y1_etf, x1_etf ) - 11.25 * piby1 / 180.0;
563 if ( fi_etf_1 < 0 ) { fi_etf_1 = pi2 + fi_etf_1; }
564
565 int Etfid_1 = (int)( fi_etf_1 / piby18 );
566 if ( Etfid_1 > 35 ) { Etfid_1 = Etfid_1 - 36; }
567
568 if ( Etfid_1 % 2 == 1 )
569 {
574 }
575 else
576 {
577 double Z_etf_2 = -136.573;
578 double Pathl_2 = Z_etf_2 /
sin( atan(
Tanl ) );
579 double phi0_2 = -( Z_etf_2 /
Tanl ) / rho;
580 double phi_2 = -( Z_etf_2 -
Dz ) /
Tanl / rho;
581 double phi1_2 = ( Z_etf_2 *
Kappa * 0.003 ) /
Tanl;
582 double phi2_2 = ( Z_etf_2 -
Dz ) *
Kappa * 0.003 /
Tanl;
583
588 int r_etf_2 = sqrt( x2_etf * x2_etf + y2_etf * y2_etf );
589
590 double x_etf_2 = helix1.x( phi_2 ).x();
591 double y_etf_2 = helix1.x( phi_2 ).y();
592 double z_etf_2 = helix1.x( phi_2 ).z();
593
594
595 double fi_etf_2 = atan2( y2_etf, x2_etf ) - 11.25 * piby1 / 180.0;
596 if ( fi_etf_2 < 0 ) { fi_etf_2 = pi2 + fi_etf_2; }
597
598 int Etfid_2 = (int)( fi_etf_2 / piby18 );
599 if ( Etfid_2 > 35 ) { Etfid_2 = Etfid_2 - 36; }
600
601 if ( Etfid_2 % 2 == 1 )
602 {
603 int tmp1 = (int)( ( fi_etf_2 + 0.5 * piby18 ) / piby18 );
604 int tmp2 = (int)( ( fi_etf_2 - 0.5 * piby18 ) / piby18 );
605 if ( tmp1 == Etfid_2 ) { Etfid_2 = tmp2; }
606 else { Etfid_2 = tmp1; }
607 }
608
613 }
614
615 return ( 2 );
616 }
617 }
618
620 {
621
622 if ( _debug )
623 {
624 printf( "\n TofFz_helix> Trk=%3d params(dr,phi,kappa,dz,tanl)="
625 "(%5.1f,%6.3f,%6.4f,%6.1f,%6.3f) R_tof %5.1f\n",
627
628 printf( " TofFz_helix> rho=%8.1f, (xc,yc)=(%8.1f,%8.1f)"
629 " (nx,ny)=(%5.2f,%5.2f)\n",
630 rho, xc, yc, nx, ny );
631
632 printf( " TofFz_helix> tof (x,y)=(%5.1f,%5.1f) and (%5.1f,%5.1f)\n", xx[0], yy[0], xx[1],
633 yy[1] );
634
635 printf( " TofFz_helix> tofid=%3d, fitof=%6.3f, w=%5.3f"
636 " (x,y,z)=(%5.1f,%5.1f,%5.1f) pathl=%5.1f cm path=%5.1f cm\n",
638 }
639 return ( -7 );
640 }
641
642 std::cerr << "TofFz_helix> Error: should not reach here" << std::endl;
643 exit( 1 );
644}
HepGeom::Point3D< double > HepPoint3D
double sin(const BesAngle a)
double cos(const BesAngle a)