BOSS 7.0.6
BESIII Offline Software System
Loading...
Searching...
No Matches
Toffz_helix.cxx
Go to the documentation of this file.
1////////////////////////////////////////////////////////////
2/// Obtains TOF hit IDs that match MDC tracks found by ///
3/// Belle Fast Tracker, fzisan. ///
4/// ///
5/// How: ///
6/// Extrapolates fzisan tracks (helix) to the inner ///
7/// surface of TOF counters to obtain corresponding ///
8/// IDs. ///
9/// ///
10/// Function: TofFz_get(double rtof, int id) ///
11/// rtof: Input TOF counter radius ///
12/// id: Input fzisan track Id ///
13/// Return IDs: ///
14/// 1: OK ///
15/// -1: wrong fzisan ID ///
16/// -3: track multiplicity <= 0 ///
17/// -5: inclomplete fitting in fzisan ///
18/// -7: bad path length or Z_tof ///
19/// ///
20/// Based on: Tof_helix class ///
21/// ///
22/// Author: H.Kichimi (for TRASAN tracks) ///
23/// ///
24/// Modifications: ///
25/// Modified for fzisan tracks S.Behari 15/01/2000 ///
26/////////////////////////////////////////////////////////////
27
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"
38
39#include "MdcGeomSvc/IMdcGeomSvc.h"
40#include "MdcGeomSvc/MdcGeoWire.h"
41#include "MdcGeomSvc/MdcGeoLayer.h"
42#include <vector>
43#include <iostream>
44#include <cstdio>
45#include <math.h>
46
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"
52
53#ifndef ENABLE_BACKWARDS_COMPATIBILITY
54 typedef HepGeom::Point3D<double> HepPoint3D;
55#endif
56
57using CLHEP::HepVector;
58using CLHEP::Hep3Vector;
59
60
61TofFz_helix::TofFz_helix(void) {
62 piby1=3.141593;
63 pi2=2.0*piby1;
64 piby44=piby1/44.0;
65 piby24=piby1/24.0;
66 piby18=piby1/18.0;
67 _debug=1;
68 _pathl_cut= 500.0;
69 _ztof_cutm=-140.0;
70 _ztof_cutx= 140.0;
71
72
73}
74
75
76int TofFz_helix::TofFz_Get(double rtof, int id, double fzisan[]) {
77 // crossing points of the helix on the tof cylinder
78 double xc,yc; // center of the circle
79 double rho; // radius of tof cyclinder
80 double xx[2],yy[2]; // coordinates of hits
81 double tofx,tofy; // (x,y) on the radius of R_tof
82 double nx,ny; // direction cosines of mom vector at vertex
83 double cosx,sinx;
84 double aa,ca,cb,cc,detm;
85
86 // fzisan track ID and TOF radius
87 int Id_cdfztrk = id; R_tof = rtof;
88
89 if( Id_cdfztrk<0 ) {
90 if (_debug) std::cout << " TofFz_helix *** wrong id_cdfztrk ="
91 << Id_cdfztrk << std::endl;
92 return(-1);
93 }
94
95 const HepPoint3D pivot1(0.,0.,0.);
96 HepVector a(5,0);
97 a[0] = fzisan[0];
98 a[1] = fzisan[1];
99 a[2] = fzisan[2];
100 a[3] = fzisan[3];
101 a[4] = fzisan[4];
102
103
104 // Ill-fitted track in fzisan (dz=tanl=-9999.)
105 if (abs(a[3])>50.0 || abs(a[4])>500.0) return (-5);
106 // if (abs(a[3])>10.0 || abs(a[4])>500.0) return (-5);
107
108 // D e f i n e h e l i x
109 Helix helix1(pivot1,a);
110 Dr = helix1.a()[0];
111 Phi0 = helix1.a()[1];
112 Kappa = helix1.a()[2];
113 Dz = helix1.a()[3];
114 Tanl = helix1.a()[4];
115
116 // check cdfztrk hit on tof cyclinder
117 rho = helix1.radius();
118 // cout<<" Toffz_helix:: rho ="<<rho<<endl; // radius of the circle in cm
119 HepPoint3D xyz = helix1.center();
120 xc = xyz(0);
121 yc = xyz(1);
122 // cout<<"Toffz_helix::helix center:xc="<<xc<<",yc= "<<yc<<endl;
123
124 Hep3Vector vxyz = helix1.direction();
125 nx = vxyz(0); // direction of track at the vertex
126 ny = vxyz(1);
127 // cout<<"Toffz_helix::direction: nx= "<<nx<<", ny="<<ny<<endl;
128
129 if(fabs(rho)>R_tof/2){
130
131 // get hit on tof cylinder at radius R_tof;
132
133 if( xc==0.0 && yc==0.0 ) return(-3);
134 // coefficients of quadratic equation
135 ca = xc*xc + yc*yc ;
136 aa = 0.5 * ( ca - rho*rho + R_tof*R_tof );
137
138 if ( xc != 0.0 ) {
139 cb = aa * yc;
140 cc = aa*aa - R_tof*R_tof*xc*xc;
141 // determinant
142 detm = cb*cb - ca*cc;
143 if ( detm > 0.0 ) {
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;
148 } else return(-1);
149 }
150 else{
151 cb = aa * xc;
152 cc = aa*aa - R_tof*R_tof*yc*yc;
153 // determinant
154 detm = cb*cb - ca*cc;
155 if ( detm > 0.0 ) {
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;
160 } else return(-2);
161 }
162
163 // choose one of hits according to track direction at the vertex.
164
165 if( xx[0]*nx + yy[0]*ny > 0.0 ) { tofx = xx[0]; tofy = yy[0]; }
166 else { tofx = xx[1]; tofy = yy[1]; }
167
168 double fi = atan2(tofy ,tofx ); // atna2 from -pi to pi
169
170
171 if( fi < 0.0 ) fi=pi2+fi;
172 Fi_tof = fi;
173 // Tofid = (int) ( Fi_tof/piby44 + 1.0 );
174 Tofid = (int) ( Fi_tof/piby44 );
175
176 // get phi and z on the tof counter
177 const HepPoint3D pivot2( tofx, tofy, 0.0);
178 helix1.pivot(pivot2);
179
180 // corrected by Ozaki san to avoid negative path length on Aug.10,99
181 Phi1=helix1.a()[1];
182 Z_tof=helix1.a()[3];
183
184 double dphi = (-xc*(tofx-xc)-yc*(tofy-yc))/sqrt(xc*xc+yc*yc)/fabs(rho);
185 dphi = acos(dphi);
186 Pathl = fabs(rho*dphi);
187 double coscor=sqrt(1.0+Tanl*Tanl);
188 Pathl=Pathl*coscor;
189
190 // path length in tof scintillator
191 Hep3Vector vxyz1 = helix1.direction();
192 nx = vxyz1(0); // direction of track at the vertex
193 ny = vxyz1(1);
194 // cout<<" Toffz_helix::track direction: nx="<<nx<<",ny="<<ny<<endl;
195 double corxy=(nx*tofx+ny*tofy)/R_tof;
196 Path_tof=4.0*coscor/corxy;
197 // cout<<" Toffz_helix::Path_tof="<< Path_tof<<endl;
198 if(abs(Z_tof)<115) return (1);
199 // if(Z_tof>116.5){
200
201
202 //std::cout << " helix1.a()[3] " << helix1.a()[3] <<std::endl;
203
204 if(Z_tof>=115){
205 // Scintillator East Endcap TOF
206 Z_tof = 133;
207 Pathl = Z_tof/sin(atan(Tanl));
208 double phi0 = -(Z_tof/Tanl)/rho;
209 double phi = -(Z_tof-Dz)/Tanl/rho;
210 double phi1 = (Z_tof*Kappa*0.003)/Tanl;
211 double phi2 = (Z_tof-Dz)*Kappa*0.003/Tanl;
212
213 double x_endtof = Dr*cos(Phi0)+(1/(Kappa*0.003))*(cos(Phi0)-cos(Phi0+phi));
214 double y_endtof = Dr*sin(Phi0)+(1/(Kappa*0.003))*(sin(Phi0)-sin(Phi0+phi));
215 r_endtof = sqrt(x_endtof * x_endtof + y_endtof * y_endtof);
216
217 Helix helix1(pivot1,a);
218
219 double x1_endtof=helix1.x(phi).x();
220 double y1_endtof=helix1.x(phi).y();
221 double z1_endtof=helix1.x(phi).z();
222
223 double fi_endtof = atan2(y_endtof,x_endtof ); // atna2 from -pi to pi
224 if (fi_endtof<0) fi_endtof=pi2+fi_endtof;
225
226 Tofid =(int)(fi_endtof/piby24);
227 if(Tofid>47)Tofid=Tofid-48;
228
229 // MRPC East Endcap TOF
230 // using hardware provide number: 1330+1+1+4.73 / mingming 133.7
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;
236 double phi2_1 = (Z_etf_1-Dz)*Kappa*0.003/Tanl;
237
238 double x1_etf = Dr*cos(Phi0)+(1/(Kappa*0.003))*(cos(Phi0)-cos(Phi0+phi_1));
239 double y1_etf = Dr*sin(Phi0)+(1/(Kappa*0.003))*(sin(Phi0)-sin(Phi0+phi_1));
240 double r_etf_1 = sqrt( x1_etf*x1_etf + y1_etf*y1_etf );
241
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();
245
246 // the sysmetric axis of east endcap has a degree shift
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; }
249
250 int Etfid_1 = (int)(fi_etf_1/piby18);
251 if( Etfid_1>35 ) { Etfid_1 = Etfid_1 - 36; }
252
253 // Inner Layer
254 if( Etfid_1%2 == 1 ) {
255 Etfid = Etfid_1;
256 r_etf = r_etf_1;
257 Z_etf = Z_etf_1;
258 Path_etf = Pathl_1;
259 }
260 // Outer Layer
261 else {
262 // using hardware provide number 1357.5+1+2.5+4.73 / mingming 136.4
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;
268 double phi2_2 = (Z_etf_2-Dz)*Kappa*0.003/Tanl;
269
270 double x2_etf = Dr*cos(Phi0)+(1/(Kappa*0.003))*(cos(Phi0)-cos(Phi0+phi_2));
271 double y2_etf = Dr*sin(Phi0)+(1/(Kappa*0.003))*(sin(Phi0)-sin(Phi0+phi_2));
272 double r_etf_2 = sqrt( x2_etf*x2_etf + y2_etf*y2_etf );
273
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();
277
278 // the sysmetric axis of east endcap has a 3.75 degree shift
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; }
281
282 int Etfid_2 = (int)(fi_etf_2/piby18);
283 if( Etfid_2>35 ) { Etfid_2 = Etfid_2 - 36; }
284
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; }
290 }
291
292 Etfid = Etfid_2;
293 r_etf = r_etf_2;
294 Z_etf = Z_etf_2;
295 Path_etf = Pathl_2;
296 }
297
298 return (0);
299
300 }
301 else if (Z_tof<=-115){
302 // Scintillator West Endcap TOF
303 Z_tof = -133;
304 Pathl = Z_tof/sin(atan(Tanl));
305 double phi0 = -(Z_tof/Tanl)/rho;
306 double phi = -(Z_tof-Dz)/Tanl/rho;
307 double phi1 = (Z_tof*Kappa*0.003)/Tanl;
308 double phi2 = (Z_tof-Dz)*Kappa*0.003/Tanl;
309
310 double x_endtof = Dr*cos(Phi0)+(1/(Kappa*0.003))*(cos(Phi0)-cos(Phi0+phi));
311 double y_endtof = Dr*sin(Phi0)+(1/(Kappa*0.003))*(sin(Phi0)-sin(Phi0+phi));
312 r_endtof = sqrt(x_endtof * x_endtof + y_endtof * y_endtof);
313
314 Helix helix1(pivot1,a);
315
316 double x1_endtof=helix1.x(phi).x();
317 double y1_endtof=helix1.x(phi).y();
318 double z1_endtof=helix1.x(phi).z();
319
320 double fi_endtof = atan2(y_endtof,x_endtof ); // atna2 from -pi to pi
321 if (fi_endtof<0) fi_endtof=pi2+fi_endtof;
322
323 Tofid =(int)(fi_endtof/piby24);
324 if(Tofid>47) Tofid=Tofid-48;
325
326 // MRPC West Endcap TOF
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;
332 double phi2_1 = (Z_etf_1-Dz)*Kappa*0.003/Tanl;
333
334 double x1_etf = Dr*cos(Phi0)+(1/(Kappa*0.003))*(cos(Phi0)-cos(Phi0+phi_1));
335 double y1_etf = Dr*sin(Phi0)+(1/(Kappa*0.003))*(sin(Phi0)-sin(Phi0+phi_1));
336 double r_etf_1 = sqrt( x1_etf*x1_etf + y1_etf*y1_etf );
337
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();
341
342 // the sysmetric axis of west endcap has a 16.25 degree shift
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; }
345
346 int Etfid_1 = (int)(fi_etf_1/piby18);
347 if( Etfid_1>35 ) { Etfid_1 = Etfid_1 - 36; }
348
349 // Inner Layer
350 if( Etfid_1%2 == 1 ) {
351 Etfid = Etfid_1;
352 r_etf = r_etf_1;
353 Z_etf = Z_etf_1;
354 Path_etf = Pathl_1;
355 }
356 // Outer Layer
357 else {
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;
363 double phi2_2 = (Z_etf_2-Dz)*Kappa*0.003/Tanl;
364
365 double x2_etf = Dr*cos(Phi0)+(1/(Kappa*0.003))*(cos(Phi0)-cos(Phi0+phi_2));
366 double y2_etf = Dr*sin(Phi0)+(1/(Kappa*0.003))*(sin(Phi0)-sin(Phi0+phi_2));
367 double r_etf_2 = sqrt( x2_etf*x2_etf + y2_etf*y2_etf );
368
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();
372
373 // the sysmetric axis of west endcap has a 16.25 degree shift
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; }
376
377 int Etfid_2 = (int)(fi_etf_2/piby18);
378 if( Etfid_2>35 ) { Etfid_2 = Etfid_2 - 36; }
379
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; }
385 }
386
387 Etfid = Etfid_2;
388 r_etf = r_etf_2;
389 Z_etf = Z_etf_2;
390 Path_etf = Pathl_2;
391 }
392
393 return (2);
394
395 }
396
397 }
398 else {
399 if(Tanl>0){
400 // Scintillator East Endcap TOF
401 Z_tof = 133;
402 Pathl = Z_tof/sin(atan(Tanl));
403 double phi0 = -(Z_tof/Tanl)/rho;
404 double phi = -(Z_tof-Dz)/Tanl/rho;
405 double phi1 = (Z_tof*Kappa*0.003)/Tanl;
406 double phi2 = (Z_tof-Dz)*Kappa*0.003/Tanl;
407
408 double x_endtof=Dr*cos(Phi0)+(1/(Kappa*0.003))*(cos(Phi0)-cos(Phi0+phi));
409 double y_endtof=Dr*sin(Phi0)+(1/(Kappa*0.003))*(sin(Phi0)-sin(Phi0+phi));
410 r_endtof=sqrt(x_endtof * x_endtof + y_endtof * y_endtof);
411
412 double fi_endtof = atan2(y_endtof,x_endtof ); // atna2 from -pi to pi
413 if (fi_endtof<0) fi_endtof=pi2+fi_endtof;
414 Tofid =(int)(fi_endtof/piby24);
415 if(Tofid>47)Tofid=Tofid-48;
416
417 // MRPC East Endcap TOF
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;
423 double phi2_1 = (Z_etf_1-Dz)*Kappa*0.003/Tanl;
424
425 double x1_etf = Dr*cos(Phi0)+(1/(Kappa*0.003))*(cos(Phi0)-cos(Phi0+phi_1));
426 double y1_etf = Dr*sin(Phi0)+(1/(Kappa*0.003))*(sin(Phi0)-sin(Phi0+phi_1));
427 double r_etf_1 = sqrt( x1_etf*x1_etf + y1_etf*y1_etf );
428
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();
432
433 // the sysmetric axis of east endcap has a 3.75 degree shift
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; }
436
437 int Etfid_1 = (int)(fi_etf_1/piby18);
438 if( Etfid_1>35 ) { Etfid_1 = Etfid_1 - 36; }
439
440 // Inner Layer
441 if( Etfid_1%2 == 1 ) {
442 Etfid = Etfid_1;
443 r_etf = r_etf_1;
444 Z_etf = Z_etf_1;
445 Path_etf = Pathl_1;
446 }
447 // Outer Layer
448 else {
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;
454 double phi2_2 = (Z_etf_2-Dz)*Kappa*0.003/Tanl;
455
456 double x2_etf = Dr*cos(Phi0)+(1/(Kappa*0.003))*(cos(Phi0)-cos(Phi0+phi_2));
457 double y2_etf = Dr*sin(Phi0)+(1/(Kappa*0.003))*(sin(Phi0)-sin(Phi0+phi_2));
458 double r_etf_2 = sqrt( x2_etf*x2_etf + y2_etf*y2_etf );
459
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();
463
464 // the sysmetric axis of east endcap has a 3.75 degree shift
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; }
467
468 int Etfid_2 = (int)(fi_etf_2/piby18);
469 if( Etfid_2>35 ) { Etfid_2 = Etfid_2 - 36; }
470
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; }
476 }
477 Etfid = Etfid_2;
478 r_etf = r_etf_2;
479 Z_etf = Z_etf_2;
480 Path_etf = Pathl_2;
481 }
482
483 return (0);
484 }
485 else{
486 // Scintillator West Endcap TOF
487 Z_tof = -133;
488 Pathl = Z_tof/sin(atan(Tanl));
489 double phi0 = -(Z_tof/Tanl)/rho;
490 double phi = -(Z_tof-Dz)/Tanl/rho;
491 double phi1 = (Z_tof*Kappa*0.003)/Tanl;
492 double phi2 = (Z_tof-Dz)*Kappa*0.003/Tanl;
493
494 double x_endtof=Dr*cos(Phi0)+(1/(Kappa*0.003))*(cos(Phi0)-cos(Phi0+phi));
495 double y_endtof=Dr*sin(Phi0)+(1/(Kappa*0.003))*(sin(Phi0)-sin(Phi0+phi));
496 r_endtof=sqrt(x_endtof * x_endtof + y_endtof * y_endtof);
497
498 double fi_endtof = atan2(y_endtof,x_endtof ); // atna2 from -pi to pi
499 if (fi_endtof<0) fi_endtof=pi2+fi_endtof;
500 Tofid =(int)(fi_endtof/piby24);
501 if(Tofid>47)Tofid=Tofid-48;
502
503 // MRPC West Endcap TOF
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;
509 double phi2_1 = (Z_etf_1-Dz)*Kappa*0.003/Tanl;
510
511 double x1_etf = Dr*cos(Phi0)+(1/(Kappa*0.003))*(cos(Phi0)-cos(Phi0+phi_1));
512 double y1_etf = Dr*sin(Phi0)+(1/(Kappa*0.003))*(sin(Phi0)-sin(Phi0+phi_1));
513 double r_etf_1 = sqrt( x1_etf*x1_etf + y1_etf*y1_etf );
514
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();
518
519 // the sysmetric axis of west endcap has a 16.25 degree shift
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; }
522
523 int Etfid_1 = (int)(fi_etf_1/piby18);
524 if( Etfid_1>35 ) { Etfid_1 = Etfid_1 - 36; }
525
526 if( Etfid_1%2 == 1 ) {
527 Etfid = Etfid_1;
528 r_etf = r_etf_1;
529 Z_etf = Z_etf_1;
530 Path_etf = Pathl_1;
531 }
532 else {
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;
538 double phi2_2 = (Z_etf_2-Dz)*Kappa*0.003/Tanl;
539
540 double x2_etf = Dr*cos(Phi0)+(1/(Kappa*0.003))*(cos(Phi0)-cos(Phi0+phi_2));
541 double y2_etf = Dr*sin(Phi0)+(1/(Kappa*0.003))*(sin(Phi0)-sin(Phi0+phi_2));
542 int r_etf_2 = sqrt( x2_etf*x2_etf + y2_etf*y2_etf );
543
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();
547
548 // the sysmetric axis of west endcap has a 16.25 degree shift
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; }
551
552 int Etfid_2 = (int)(fi_etf_2/piby18);
553 if( Etfid_2>35 ) { Etfid_2 = Etfid_2 - 36; }
554
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; }
560 }
561
562 Etfid = Etfid_2;
563 r_etf = r_etf_2;
564 Z_etf = Z_etf_2;
565 Path_etf = Pathl_2;
566 }
567
568 return (2);
569 }
570 }
571
572 if (abs(Pathl) > _pathl_cut ||
573 Z_tof < _ztof_cutm || Z_tof > _ztof_cutx) {
574 // Bad path length or Z_tof
575 if (_debug) {
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",
578 Id_cdfztrk,Dr,Phi0,Kappa,Dz,Tanl,R_tof);
579
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);
582
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]);
585
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",
588 Tofid,Fi_tof,W_tof,tofx,tofy,Z_tof,Pathl,Path_tof);
589 }
590 return (-7);
591 }
592
593}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: Toffz_helix.cxx:54
Toffz_helix.h
Helix::center
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/TrackUtil/Helix.h:203
Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/src/Helix.cxx:129
Helix::radius
double radius(void) const
returns radious of helix.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/TrackUtil/Helix.h:215
Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/TrackUtil/Helix.h:263
Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/TrackUtil/Helix.h:209
Helix::direction
Hep3Vector direction(double dPhi=0.) const
returns direction vector after rotating angle dPhi in phi direction.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/TrackUtil/Helix.h:221
TofFz_helix::TofFz_Get
int TofFz_Get(double, int, double[])
Definition: Toffz_helix.cxx:76
TofFz_helix::Kappa
double Kappa
Definition: Toffz_helix.h:39
TofFz_helix::Path_tof
double Path_tof
Definition: Toffz_helix.h:33
TofFz_helix::r_endtof
double r_endtof
Definition: Toffz_helix.h:42
TofFz_helix::Pathl
double Pathl
Definition: Toffz_helix.h:32
TofFz_helix::Path_etf
double Path_etf
Definition: Toffz_helix.h:34
TofFz_helix::TofFz_helix
TofFz_helix(void)
Definition: Toffz_helix.cxx:61
TofFz_helix::Dr
double Dr
Definition: Toffz_helix.h:39
TofFz_helix::r_etf
double r_etf
Definition: Toffz_helix.h:43
TofFz_helix::Dz
double Dz
Definition: Toffz_helix.h:39
TofFz_helix::Phi0
double Phi0
Definition: Toffz_helix.h:39
TofFz_helix::Tofid
int Tofid
Definition: Toffz_helix.h:27
TofFz_helix::W_tof
double W_tof
Definition: Toffz_helix.h:31
TofFz_helix::Etfid
int Etfid
Definition: Toffz_helix.h:28
TofFz_helix::Tanl
double Tanl
Definition: Toffz_helix.h:39
TofFz_helix::R_tof
double R_tof
Definition: Toffz_helix.h:29
TofFz_helix::Phi1
double Phi1
Definition: Toffz_helix.h:41
TofFz_helix::Z_tof
double Z_tof
Definition: Toffz_helix.h:35
TofFz_helix::Z_etf
double Z_etf
Definition: Toffz_helix.h:36
TofFz_helix::Fi_tof
double Fi_tof
Definition: Toffz_helix.h:30