28#include "GaudiKernel/MsgStream.h"
29#include "GaudiKernel/AlgFactory.h"
30#include "GaudiKernel/ISvcLocator.h"
31#include "GaudiKernel/SmartDataPtr.h"
32#include "GaudiKernel/IDataProviderSvc.h"
33#include "GaudiKernel/PropertyMgr.h"
34#include "EventModel/Event.h"
35#include "MdcRawEvent/MdcDigi.h"
36#include "TofRawEvent/TofDigi.h"
37#include "McTruth/McParticle.h"
39#include "MdcGeomSvc/IMdcGeomSvc.h"
40#include "MdcGeomSvc/MdcGeoWire.h"
41#include "MdcGeomSvc/MdcGeoLayer.h"
47#include "EsTimeAlg/Toffz_helix.h"
48#include "CLHEP/Vector/ThreeVector.h"
49#include "MdcFastTrkAlg/FTFinder.h"
50#include "TrackUtil/Helix.h"
51#include "CLHEP/Geometry/Point3D.h"
53#ifndef ENABLE_BACKWARDS_COMPATIBILITY
57using CLHEP::HepVector;
58using CLHEP::Hep3Vector;
84 double aa,ca,cb,cc,detm;
87 int Id_cdfztrk = id;
R_tof = rtof;
90 if (_debug) std::cout <<
" TofFz_helix *** wrong id_cdfztrk ="
91 << Id_cdfztrk << std::endl;
105 if (
abs(a[3])>50.0 ||
abs(a[4])>500.0)
return (-5);
109 Helix helix1(pivot1,a);
111 Phi0 = helix1.
a()[1];
114 Tanl = helix1.
a()[4];
129 if(fabs(rho)>
R_tof/2){
133 if( xc==0.0 && yc==0.0 )
return(-3);
142 detm = cb*cb - ca*cc;
144 yy[0] = ( cb + sqrt(detm) )/ ca;
145 xx[0] = ( aa - yy[0]*yc )/xc;
146 yy[1] = ( cb - sqrt(detm) )/ ca;
147 xx[1] = ( aa - yy[1]*yc )/xc;
154 detm = cb*cb - ca*cc;
156 xx[0] = ( cb + sqrt(detm) )/ca;
157 yy[0] = ( aa - xx[0]*xc )/yc;
158 xx[1] = ( cb - sqrt(detm) )/ca;
159 yy[1] = ( aa - xx[1]*xc )/yc;
165 if( xx[0]*nx + yy[0]*ny > 0.0 ) { tofx = xx[0]; tofy = yy[0]; }
166 else { tofx = xx[1]; tofy = yy[1]; }
168 double fi = atan2(tofy ,tofx );
171 if( fi < 0.0 ) fi=pi2+fi;
178 helix1.
pivot(pivot2);
184 double dphi = (-xc*(tofx-xc)-yc*(tofy-yc))/sqrt(xc*xc+yc*yc)/fabs(rho);
186 Pathl = fabs(rho*dphi);
195 double corxy=(nx*tofx+ny*tofy)/
R_tof;
215 r_endtof = sqrt(x_endtof * x_endtof + y_endtof * y_endtof);
217 Helix helix1(pivot1,a);
219 double x1_endtof=helix1.
x(phi).x();
220 double y1_endtof=helix1.
x(phi).y();
221 double z1_endtof=helix1.
x(phi).z();
223 double fi_endtof = atan2(y_endtof,x_endtof );
224 if (fi_endtof<0) fi_endtof=pi2+fi_endtof;
226 Tofid =(int)(fi_endtof/piby24);
231 double Z_etf_1 = 133.673;
232 double Pathl_1 = Z_etf_1/
sin(atan(
Tanl));
233 double phi0_1 = -(Z_etf_1/
Tanl)/rho;
234 double phi_1 = -(Z_etf_1-
Dz)/
Tanl/rho;
235 double phi1_1 = (Z_etf_1*
Kappa*0.003)/
Tanl;
240 double r_etf_1 = sqrt( x1_etf*x1_etf + y1_etf*y1_etf );
242 double x_etf_1 = helix1.
x(phi_1).x();
243 double y_etf_1 = helix1.
x(phi_1).y();
244 double z_etf_1 = helix1.
x(phi_1).z();
247 double fi_etf_1 = atan2( y1_etf, x1_etf ) + 1.25*piby1/180.0;
248 if( fi_etf_1<0 ) { fi_etf_1 = pi2 + fi_etf_1; }
250 int Etfid_1 = (int)(fi_etf_1/piby18);
251 if( Etfid_1>35 ) { Etfid_1 = Etfid_1 - 36; }
254 if( Etfid_1%2 == 1 ) {
263 double Z_etf_2 = 136.573;
264 double Pathl_2 = Z_etf_2/
sin(atan(
Tanl));
265 double phi0_2 = -(Z_etf_2/
Tanl)/rho;
266 double phi_2 = -(Z_etf_2-
Dz)/
Tanl/rho;
267 double phi1_2 = (Z_etf_2*
Kappa*0.003)/
Tanl;
272 double r_etf_2 = sqrt( x2_etf*x2_etf + y2_etf*y2_etf );
274 double x_etf_2 = helix1.
x(phi_2).x();
275 double y_etf_2 = helix1.
x(phi_2).y();
276 double z_etf_2 = helix1.
x(phi_2).z();
279 double fi_etf_2 = atan2( y2_etf, x2_etf ) + 1.25*piby1/180.0 ;
280 if( fi_etf_2<0 ) { fi_etf_2 = pi2 + fi_etf_2; }
282 int Etfid_2 = (int)(fi_etf_2/piby18);
283 if( Etfid_2>35 ) { Etfid_2 = Etfid_2 - 36; }
285 if( Etfid_2%2 == 1 ) {
286 int tmp1 = (int)((fi_etf_2+0.5*piby18)/piby18);
287 int tmp2 = (int)((fi_etf_2-0.5*piby18)/piby18);
288 if( tmp1 == Etfid_2 ) { Etfid_2 = tmp2; }
289 else { Etfid_2 = tmp1; }
301 else if (
Z_tof<=-115){
312 r_endtof = sqrt(x_endtof * x_endtof + y_endtof * y_endtof);
314 Helix helix1(pivot1,a);
316 double x1_endtof=helix1.
x(phi).x();
317 double y1_endtof=helix1.
x(phi).y();
318 double z1_endtof=helix1.
x(phi).z();
320 double fi_endtof = atan2(y_endtof,x_endtof );
321 if (fi_endtof<0) fi_endtof=pi2+fi_endtof;
323 Tofid =(int)(fi_endtof/piby24);
327 double Z_etf_1 = -133.673;
328 double Pathl_1 = Z_etf_1/
sin(atan(
Tanl));
329 double phi0_1 = -(Z_etf_1/
Tanl)/rho;
330 double phi_1 = -(Z_etf_1-
Dz)/
Tanl/rho;
331 double phi1_1 = (Z_etf_1*
Kappa*0.003)/
Tanl;
336 double r_etf_1 = sqrt( x1_etf*x1_etf + y1_etf*y1_etf );
338 double x_etf_1 = helix1.
x(phi_1).x();
339 double y_etf_1 = helix1.
x(phi_1).y();
340 double z_etf_1 = helix1.
x(phi_1).z();
343 double fi_etf_1 = atan2( y1_etf, x1_etf ) - 11.25*piby1/180.0;
344 if( fi_etf_1<0 ) { fi_etf_1 = pi2 + fi_etf_1; }
346 int Etfid_1 = (int)(fi_etf_1/piby18);
347 if( Etfid_1>35 ) { Etfid_1 = Etfid_1 - 36; }
350 if( Etfid_1%2 == 1 ) {
358 double Z_etf_2 = -136.573;
359 double Pathl_2 = Z_etf_2/
sin(atan(
Tanl));
360 double phi0_2 = -(Z_etf_2/
Tanl)/rho;
361 double phi_2 = -(Z_etf_2-
Dz)/
Tanl/rho;
362 double phi1_2 = (Z_etf_2*
Kappa*0.003)/
Tanl;
367 double r_etf_2 = sqrt( x2_etf*x2_etf + y2_etf*y2_etf );
369 double x_etf_2 = helix1.
x(phi_2).x();
370 double y_etf_2 = helix1.
x(phi_2).y();
371 double z_etf_2 = helix1.
x(phi_2).z();
374 double fi_etf_2 = atan2( y2_etf, x2_etf ) - 11.25*piby1/180.0;
375 if( fi_etf_2<0 ) { fi_etf_2 = pi2 + fi_etf_2; }
377 int Etfid_2 = (int)(fi_etf_2/piby18);
378 if( Etfid_2>35 ) { Etfid_2 = Etfid_2 - 36; }
380 if( Etfid_2%2 == 1 ) {
381 int tmp1 = (int)((fi_etf_2+0.5*piby18)/piby18);
382 int tmp2 = (int)((fi_etf_2-0.5*piby18)/piby18);
383 if( tmp1 == Etfid_2 ) { Etfid_2 = tmp2; }
384 else { Etfid_2 = tmp1; }
410 r_endtof=sqrt(x_endtof * x_endtof + y_endtof * y_endtof);
412 double fi_endtof = atan2(y_endtof,x_endtof );
413 if (fi_endtof<0) fi_endtof=pi2+fi_endtof;
414 Tofid =(int)(fi_endtof/piby24);
418 double Z_etf_1 = 133.673;
419 double Pathl_1 = Z_etf_1/
sin(atan(
Tanl));
420 double phi0_1 = -(Z_etf_1/
Tanl)/rho;
421 double phi_1 = -(Z_etf_1-
Dz)/
Tanl/rho;
422 double phi1_1 = (Z_etf_1*
Kappa*0.003)/
Tanl;
427 double r_etf_1 = sqrt( x1_etf*x1_etf + y1_etf*y1_etf );
429 double x_etf_1 = helix1.
x(phi_1).x();
430 double y_etf_1 = helix1.
x(phi_1).y();
431 double z_etf_1 = helix1.
x(phi_1).z();
434 double fi_etf_1 = atan2( y1_etf, x1_etf ) + 1.25*piby1/180.0;
435 if( fi_etf_1<0 ) { fi_etf_1 = pi2 + fi_etf_1; }
437 int Etfid_1 = (int)(fi_etf_1/piby18);
438 if( Etfid_1>35 ) { Etfid_1 = Etfid_1 - 36; }
441 if( Etfid_1%2 == 1 ) {
449 double Z_etf_2 = 136.573;
450 double Pathl_2 = Z_etf_2/
sin(atan(
Tanl));
451 double phi0_2 = -(Z_etf_2/
Tanl)/rho;
452 double phi_2 = -(Z_etf_2-
Dz)/
Tanl/rho;
453 double phi1_2 = (Z_etf_2*
Kappa*0.003)/
Tanl;
458 double r_etf_2 = sqrt( x2_etf*x2_etf + y2_etf*y2_etf );
460 double x_etf_2 = helix1.
x(phi_2).x();
461 double y_etf_2 = helix1.
x(phi_2).y();
462 double z_etf_2 = helix1.
x(phi_2).z();
465 double fi_etf_2 = atan2( y2_etf, x2_etf ) + 1.25*piby1/180.0;
466 if( fi_etf_2<0 ) { fi_etf_2 = pi2 + fi_etf_2; }
468 int Etfid_2 = (int)(fi_etf_2/piby18);
469 if( Etfid_2>35 ) { Etfid_2 = Etfid_2 - 36; }
471 if( Etfid_2%2 == 1 ) {
472 int tmp1 = (int)((fi_etf_2+0.5*piby18)/piby18);
473 int tmp2 = (int)((fi_etf_2-0.5*piby18)/piby18);
474 if( tmp1 == Etfid_2 ) { Etfid_2 = tmp2; }
475 else { Etfid_2 = tmp1; }
496 r_endtof=sqrt(x_endtof * x_endtof + y_endtof * y_endtof);
498 double fi_endtof = atan2(y_endtof,x_endtof );
499 if (fi_endtof<0) fi_endtof=pi2+fi_endtof;
500 Tofid =(int)(fi_endtof/piby24);
504 double Z_etf_1 = -133.673;
505 double Pathl_1 = Z_etf_1/
sin(atan(
Tanl));
506 double phi0_1 = -(Z_etf_1/
Tanl)/rho;
507 double phi_1 = -(Z_etf_1-
Dz)/
Tanl/rho;
508 double phi1_1 = (Z_etf_1*
Kappa*0.003)/
Tanl;
513 double r_etf_1 = sqrt( x1_etf*x1_etf + y1_etf*y1_etf );
515 double x_etf_1 = helix1.
x(phi_1).x();
516 double y_etf_1 = helix1.
x(phi_1).y();
517 double z_etf_1 = helix1.
x(phi_1).z();
520 double fi_etf_1 = atan2( y1_etf, x1_etf ) - 11.25*piby1/180.0;
521 if( fi_etf_1<0 ) { fi_etf_1 = pi2 + fi_etf_1; }
523 int Etfid_1 = (int)(fi_etf_1/piby18);
524 if( Etfid_1>35 ) { Etfid_1 = Etfid_1 - 36; }
526 if( Etfid_1%2 == 1 ) {
533 double Z_etf_2 = -136.573;
534 double Pathl_2 = Z_etf_2/
sin(atan(
Tanl));
535 double phi0_2 = -(Z_etf_2/
Tanl)/rho;
536 double phi_2 = -(Z_etf_2-
Dz)/
Tanl/rho;
537 double phi1_2 = (Z_etf_2*
Kappa*0.003)/
Tanl;
542 int r_etf_2 = sqrt( x2_etf*x2_etf + y2_etf*y2_etf );
544 double x_etf_2 = helix1.
x(phi_2).x();
545 double y_etf_2 = helix1.
x(phi_2).y();
546 double z_etf_2 = helix1.
x(phi_2).z();
549 double fi_etf_2 = atan2( y2_etf, x2_etf ) - 11.25*piby1/180.0;
550 if( fi_etf_2<0 ) { fi_etf_2 = pi2 + fi_etf_2; }
552 int Etfid_2 = (int)(fi_etf_2/piby18);
553 if( Etfid_2>35 ) { Etfid_2 = Etfid_2 - 36; }
555 if( Etfid_2%2 == 1 ) {
556 int tmp1 = (int)((fi_etf_2+0.5*piby18)/piby18);
557 int tmp2 = (int)((fi_etf_2-0.5*piby18)/piby18);
558 if( tmp1 == Etfid_2 ) { Etfid_2 = tmp2; }
559 else { Etfid_2 = tmp1; }
573 Z_tof < _ztof_cutm || Z_tof > _ztof_cutx) {
576 printf(
"\n TofFz_helix> Trk=%3d params(dr,phi,kappa,dz,tanl)="
577 "(%5.1f,%6.3f,%6.4f,%6.1f,%6.3f) R_tof %5.1f\n",
580 printf(
" TofFz_helix> rho=%8.1f, (xc,yc)=(%8.1f,%8.1f)"
581 " (nx,ny)=(%5.2f,%5.2f)\n",rho,xc,yc,nx,ny);
583 printf(
" TofFz_helix> tof (x,y)=(%5.1f,%5.1f) and (%5.1f,%5.1f)\n",
584 xx[0],yy[0],xx[1],yy[1]);
586 printf(
" TofFz_helix> tofid=%3d, fitof=%6.3f, w=%5.3f"
587 " (x,y,z)=(%5.1f,%5.1f,%5.1f) pathl=%5.1f cm path=%5.1f cm\n",
double sin(const BesAngle a)
double cos(const BesAngle a)
HepGeom::Point3D< double > HepPoint3D
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
double radius(void) const
returns radious of helix.
const HepVector & a(void) const
returns helix parameters.
const HepPoint3D & pivot(void) const
returns pivot position.
Hep3Vector direction(double dPhi=0.) const
returns direction vector after rotating angle dPhi in phi direction.
int TofFz_Get(double, int, double[])