BOSS 8.0.0
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 "CLHEP/Geometry/Point3D.h"
29#include "CLHEP/Vector/ThreeVector.h"
30#include "GaudiKernel/IDataProviderSvc.h"
31#include "GaudiKernel/ISvcLocator.h"
32#include "GaudiKernel/MsgStream.h"
33#include "GaudiKernel/PropertyMgr.h"
34#include "GaudiKernel/SmartDataPtr.h"
35#include <cstdio>
36#include <iostream>
37#include <math.h>
38#include <vector>
39
40#include "EventModel/Event.h"
41#include "McTruth/McParticle.h"
42#include "MdcGeomSvc/IMdcGeomSvc.h"
43#include "MdcGeomSvc/MdcGeoLayer.h"
44#include "MdcGeomSvc/MdcGeoWire.h"
45#include "MdcRawEvent/MdcDigi.h"
46#include "TofRawEvent/TofDigi.h"
47#include "TrackUtil/Helix.h"
48
49#include "Toffz_helix.h"
50
51#ifndef ENABLE_BACKWARDS_COMPATIBILITY
52typedef HepGeom::Point3D<double> HepPoint3D;
53#endif
54
55using CLHEP::Hep3Vector;
56using CLHEP::HepVector;
57
58TofFz_helix::TofFz_helix( void ) {
59 piby1 = 3.141593;
60 pi2 = 2.0 * piby1;
61 piby44 = piby1 / 44.0;
62 piby24 = piby1 / 24.0;
63 piby18 = piby1 / 18.0;
64 _debug = 1;
65 _pathl_cut = 500.0;
66 _ztof_cutm = -140.0;
67 _ztof_cutx = 140.0;
68}
69
70int TofFz_helix::TofFz_Get( double rtof, int id, double fzisan[] ) {
71 // crossing points of the helix on the tof cylinder
72 double xc, yc; // center of the circle
73 double rho; // radius of tof cyclinder
74 double xx[2], yy[2]; // coordinates of hits
75 double tofx, tofy; // (x,y) on the radius of R_tof
76 double nx, ny; // direction cosines of mom vector at vertex
77 double cosx, sinx;
78 double aa, ca, cb, cc, detm;
79
80 // fzisan track ID and TOF radius
81 int Id_cdfztrk = id;
82 R_tof = rtof;
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
91 const HepPoint3D pivot1( 0., 0., 0. );
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 // Ill-fitted track in fzisan (dz=tanl=-9999.)
100 if ( abs( a[3] ) > 50.0 || abs( a[4] ) > 500.0 ) return ( -5 );
101 // if (abs(a[3])>10.0 || abs(a[4])>500.0) return (-5);
102
103 // D e f i n e h e l i x
104 Helix helix1( pivot1, a );
105 Dr = helix1.a()[0];
106 Phi0 = helix1.a()[1];
107 Kappa = helix1.a()[2];
108 Dz = helix1.a()[3];
109 Tanl = helix1.a()[4];
110
111 // check cdfztrk hit on tof cyclinder
112 rho = helix1.radius();
113 // cout<<" Toffz_helix:: rho ="<<rho<<endl; // radius of the circle in cm
114 HepPoint3D xyz = helix1.center();
115 xc = xyz( 0 );
116 yc = xyz( 1 );
117 // cout<<"Toffz_helix::helix center:xc="<<xc<<",yc= "<<yc<<endl;
118
119 Hep3Vector vxyz = helix1.direction();
120 nx = vxyz( 0 ); // direction of track at the vertex
121 ny = vxyz( 1 );
122 // cout<<"Toffz_helix::direction: nx= "<<nx<<", ny="<<ny<<endl;
123
124 if ( fabs( rho ) > R_tof / 2 )
125 {
126
127 // get hit on tof cylinder at radius R_tof;
128
129 if ( xc == 0.0 && yc == 0.0 ) return ( -3 );
130 // coefficients of quadratic equation
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;
137 cc = aa * aa - R_tof * R_tof * xc * xc;
138 // determinant
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;
152 cc = aa * aa - R_tof * R_tof * yc * yc;
153 // determinant
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 // choose one of hits according to track direction at the vertex.
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 ); // atna2 from -pi to pi
179
180 if ( fi < 0.0 ) fi = pi2 + fi;
181 Fi_tof = fi;
182 // Tofid = (int) ( Fi_tof/piby44 + 1.0 );
183 Tofid = (int)( Fi_tof / piby44 );
184
185 // get phi and z on the tof counter
186 const HepPoint3D pivot2( tofx, tofy, 0.0 );
187 helix1.pivot( pivot2 );
188
189 // corrected by Ozaki san to avoid negative path length on Aug.10,99
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 );
198 Pathl = Pathl * coscor;
199
200 // path length in tof scintillator
201 Hep3Vector vxyz1 = helix1.direction();
202 nx = vxyz1( 0 ); // direction of track at the vertex
203 ny = vxyz1( 1 );
204 // cout<<" Toffz_helix::track direction: nx="<<nx<<",ny="<<ny<<endl;
205 double corxy = ( nx * tofx + ny * tofy ) / R_tof;
206 Path_tof = 4.0 * coscor / corxy;
207 // cout<<" Toffz_helix::Path_tof="<< Path_tof<<endl;
208 if ( abs( Z_tof ) < 115 ) return ( 1 );
209 // if(Z_tof>116.5){
210
211 // std::cout << " helix1.a()[3] " << helix1.a()[3] <<std::endl;
212
213 if ( Z_tof >= 115 )
214 {
215 // Scintillator East Endcap TOF
216 Z_tof = 133;
217 Pathl = Z_tof / sin( atan( Tanl ) );
218 double phi0 = -( Z_tof / Tanl ) / rho;
219 double phi = -( Z_tof - Dz ) / Tanl / rho;
220 double phi1 = ( Z_tof * Kappa * 0.003 ) / Tanl;
221 double phi2 = ( Z_tof - Dz ) * Kappa * 0.003 / Tanl;
222
223 double x_endtof =
224 Dr * cos( Phi0 ) + ( 1 / ( Kappa * 0.003 ) ) * ( cos( Phi0 ) - cos( Phi0 + phi ) );
225 double y_endtof =
226 Dr * sin( Phi0 ) + ( 1 / ( Kappa * 0.003 ) ) * ( sin( Phi0 ) - sin( Phi0 + phi ) );
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 ); // atna2 from -pi to pi
236 if ( fi_endtof < 0 ) fi_endtof = pi2 + fi_endtof;
237
238 Tofid = (int)( fi_endtof / piby24 );
239 if ( Tofid > 47 ) Tofid = Tofid - 48;
240
241 // MRPC East Endcap TOF
242 // using hardware provide number: 1330+1+1+4.73 / mingming 133.7
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 =
251 Dr * cos( Phi0 ) + ( 1 / ( Kappa * 0.003 ) ) * ( cos( Phi0 ) - cos( Phi0 + phi_1 ) );
252 double y1_etf =
253 Dr * sin( Phi0 ) + ( 1 / ( Kappa * 0.003 ) ) * ( sin( Phi0 ) - sin( Phi0 + phi_1 ) );
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 // the sysmetric axis of east endcap has a degree shift
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 // Inner Layer
268 if ( Etfid_1 % 2 == 1 )
269 {
270 Etfid = Etfid_1;
271 r_etf = r_etf_1;
272 Z_etf = Z_etf_1;
273 Path_etf = Pathl_1;
274 }
275 // Outer Layer
276 else
277 {
278 // using hardware provide number 1357.5+1+2.5+4.73 / mingming 136.4
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
286 double x2_etf = Dr * cos( Phi0 ) +
287 ( 1 / ( Kappa * 0.003 ) ) * ( cos( Phi0 ) - cos( Phi0 + phi_2 ) );
288 double y2_etf = Dr * sin( Phi0 ) +
289 ( 1 / ( Kappa * 0.003 ) ) * ( sin( Phi0 ) - sin( Phi0 + phi_2 ) );
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 // the sysmetric axis of east endcap has a 3.75 degree shift
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
311 Etfid = Etfid_2;
312 r_etf = r_etf_2;
313 Z_etf = Z_etf_2;
314 Path_etf = Pathl_2;
315 }
316
317 return ( 0 );
318 }
319 else if ( Z_tof <= -115 )
320 {
321 // Scintillator West Endcap TOF
322 Z_tof = -133;
323 Pathl = Z_tof / sin( atan( Tanl ) );
324 double phi0 = -( Z_tof / Tanl ) / rho;
325 double phi = -( Z_tof - Dz ) / Tanl / rho;
326 double phi1 = ( Z_tof * Kappa * 0.003 ) / Tanl;
327 double phi2 = ( Z_tof - Dz ) * Kappa * 0.003 / Tanl;
328
329 double x_endtof =
330 Dr * cos( Phi0 ) + ( 1 / ( Kappa * 0.003 ) ) * ( cos( Phi0 ) - cos( Phi0 + phi ) );
331 double y_endtof =
332 Dr * sin( Phi0 ) + ( 1 / ( Kappa * 0.003 ) ) * ( sin( Phi0 ) - sin( Phi0 + phi ) );
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 ); // atna2 from -pi to pi
342 if ( fi_endtof < 0 ) fi_endtof = pi2 + fi_endtof;
343
344 Tofid = (int)( fi_endtof / piby24 );
345 if ( Tofid > 47 ) Tofid = Tofid - 48;
346
347 // MRPC West Endcap TOF
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 =
356 Dr * cos( Phi0 ) + ( 1 / ( Kappa * 0.003 ) ) * ( cos( Phi0 ) - cos( Phi0 + phi_1 ) );
357 double y1_etf =
358 Dr * sin( Phi0 ) + ( 1 / ( Kappa * 0.003 ) ) * ( sin( Phi0 ) - sin( Phi0 + phi_1 ) );
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 // the sysmetric axis of west endcap has a 16.25 degree shift
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 // Inner Layer
373 if ( Etfid_1 % 2 == 1 )
374 {
375 Etfid = Etfid_1;
376 r_etf = r_etf_1;
377 Z_etf = Z_etf_1;
378 Path_etf = Pathl_1;
379 }
380 // Outer Layer
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
390 double x2_etf = Dr * cos( Phi0 ) +
391 ( 1 / ( Kappa * 0.003 ) ) * ( cos( Phi0 ) - cos( Phi0 + phi_2 ) );
392 double y2_etf = Dr * sin( Phi0 ) +
393 ( 1 / ( Kappa * 0.003 ) ) * ( sin( Phi0 ) - sin( Phi0 + phi_2 ) );
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 // the sysmetric axis of west endcap has a 16.25 degree shift
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
415 Etfid = Etfid_2;
416 r_etf = r_etf_2;
417 Z_etf = Z_etf_2;
418 Path_etf = Pathl_2;
419 }
420
421 return ( 2 );
422 }
423 }
424 else
425 {
426 if ( Tanl > 0 )
427 {
428 // Scintillator East Endcap TOF
429 Z_tof = 133;
430 Pathl = Z_tof / sin( atan( Tanl ) );
431 double phi0 = -( Z_tof / Tanl ) / rho;
432 double phi = -( Z_tof - Dz ) / Tanl / rho;
433 double phi1 = ( Z_tof * Kappa * 0.003 ) / Tanl;
434 double phi2 = ( Z_tof - Dz ) * Kappa * 0.003 / Tanl;
435
436 double x_endtof =
437 Dr * cos( Phi0 ) + ( 1 / ( Kappa * 0.003 ) ) * ( cos( Phi0 ) - cos( Phi0 + phi ) );
438 double y_endtof =
439 Dr * sin( Phi0 ) + ( 1 / ( Kappa * 0.003 ) ) * ( sin( Phi0 ) - sin( Phi0 + phi ) );
440 r_endtof = sqrt( x_endtof * x_endtof + y_endtof * y_endtof );
441
442 double fi_endtof = atan2( y_endtof, x_endtof ); // atna2 from -pi to pi
443 if ( fi_endtof < 0 ) fi_endtof = pi2 + fi_endtof;
444 Tofid = (int)( fi_endtof / piby24 );
445 if ( Tofid > 47 ) Tofid = Tofid - 48;
446
447 // MRPC East Endcap TOF
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 =
456 Dr * cos( Phi0 ) + ( 1 / ( Kappa * 0.003 ) ) * ( cos( Phi0 ) - cos( Phi0 + phi_1 ) );
457 double y1_etf =
458 Dr * sin( Phi0 ) + ( 1 / ( Kappa * 0.003 ) ) * ( sin( Phi0 ) - sin( Phi0 + phi_1 ) );
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 // the sysmetric axis of east endcap has a 3.75 degree shift
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 // Inner Layer
473 if ( Etfid_1 % 2 == 1 )
474 {
475 Etfid = Etfid_1;
476 r_etf = r_etf_1;
477 Z_etf = Z_etf_1;
478 Path_etf = Pathl_1;
479 }
480 // Outer Layer
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
490 double x2_etf = Dr * cos( Phi0 ) +
491 ( 1 / ( Kappa * 0.003 ) ) * ( cos( Phi0 ) - cos( Phi0 + phi_2 ) );
492 double y2_etf = Dr * sin( Phi0 ) +
493 ( 1 / ( Kappa * 0.003 ) ) * ( sin( Phi0 ) - sin( Phi0 + phi_2 ) );
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 // the sysmetric axis of east endcap has a 3.75 degree shift
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 }
514 Etfid = Etfid_2;
515 r_etf = r_etf_2;
516 Z_etf = Z_etf_2;
517 Path_etf = Pathl_2;
518 }
519
520 return ( 0 );
521 }
522 else
523 {
524 // Scintillator West Endcap TOF
525 Z_tof = -133;
526 Pathl = Z_tof / sin( atan( Tanl ) );
527 double phi0 = -( Z_tof / Tanl ) / rho;
528 double phi = -( Z_tof - Dz ) / Tanl / rho;
529 double phi1 = ( Z_tof * Kappa * 0.003 ) / Tanl;
530 double phi2 = ( Z_tof - Dz ) * Kappa * 0.003 / Tanl;
531
532 double x_endtof =
533 Dr * cos( Phi0 ) + ( 1 / ( Kappa * 0.003 ) ) * ( cos( Phi0 ) - cos( Phi0 + phi ) );
534 double y_endtof =
535 Dr * sin( Phi0 ) + ( 1 / ( Kappa * 0.003 ) ) * ( sin( Phi0 ) - sin( Phi0 + phi ) );
536 r_endtof = sqrt( x_endtof * x_endtof + y_endtof * y_endtof );
537
538 double fi_endtof = atan2( y_endtof, x_endtof ); // atna2 from -pi to pi
539 if ( fi_endtof < 0 ) fi_endtof = pi2 + fi_endtof;
540 Tofid = (int)( fi_endtof / piby24 );
541 if ( Tofid > 47 ) Tofid = Tofid - 48;
542
543 // MRPC West Endcap TOF
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 =
552 Dr * cos( Phi0 ) + ( 1 / ( Kappa * 0.003 ) ) * ( cos( Phi0 ) - cos( Phi0 + phi_1 ) );
553 double y1_etf =
554 Dr * sin( Phi0 ) + ( 1 / ( Kappa * 0.003 ) ) * ( sin( Phi0 ) - sin( Phi0 + phi_1 ) );
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 // the sysmetric axis of west endcap has a 16.25 degree shift
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 {
570 Etfid = Etfid_1;
571 r_etf = r_etf_1;
572 Z_etf = Z_etf_1;
573 Path_etf = Pathl_1;
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
584 double x2_etf = Dr * cos( Phi0 ) +
585 ( 1 / ( Kappa * 0.003 ) ) * ( cos( Phi0 ) - cos( Phi0 + phi_2 ) );
586 double y2_etf = Dr * sin( Phi0 ) +
587 ( 1 / ( Kappa * 0.003 ) ) * ( sin( Phi0 ) - sin( Phi0 + phi_2 ) );
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 // the sysmetric axis of west endcap has a 16.25 degree shift
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
609 Etfid = Etfid_2;
610 r_etf = r_etf_2;
611 Z_etf = Z_etf_2;
612 Path_etf = Pathl_2;
613 }
614
615 return ( 2 );
616 }
617 }
618
619 if ( abs( Pathl ) > _pathl_cut || Z_tof < _ztof_cutm || Z_tof > _ztof_cutx )
620 {
621 // Bad path length or Z_tof
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",
626 Id_cdfztrk, Dr, Phi0, Kappa, Dz, Tanl, R_tof );
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",
637 Tofid, Fi_tof, W_tof, tofx, tofy, Z_tof, Pathl, Path_tof );
638 }
639 return ( -7 );
640 }
641
642 std::cerr << "TofFz_helix> Error: should not reach here" << std::endl;
643 exit( 1 );
644}
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
phi2
Double_t phi2
Definition DataBase/tau_mode.c:8
phi1
Double_t phi1
Definition DataBase/tau_mode.c:8
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
Toffz_helix.h
Helix::center
const HepPoint3D & center(void) const
returns position of helix center(z = 0.);
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:197
Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition Reconstruction/TrackUtil/src/Helix.cxx:120
Helix::radius
double radius(void) const
returns radious of helix.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:201
Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:217
Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:199
Helix::direction
Hep3Vector direction(double dPhi=0.) const
returns direction vector after rotating angle dPhi in phi direction.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TrackUtil/Helix.h:203
TofFz_helix::TofFz_Get
int TofFz_Get(double, int, double[])
Definition Toffz_helix.cxx:70
TofFz_helix::Kappa
double Kappa
Definition Toffz_helix.h:38
TofFz_helix::Path_tof
double Path_tof
Definition Toffz_helix.h:32
TofFz_helix::r_endtof
double r_endtof
Definition Toffz_helix.h:41
TofFz_helix::Pathl
double Pathl
Definition Toffz_helix.h:31
TofFz_helix::Path_etf
double Path_etf
Definition Toffz_helix.h:33
TofFz_helix::TofFz_helix
TofFz_helix(void)
Definition Toffz_helix.cxx:58
TofFz_helix::Dr
double Dr
Definition Toffz_helix.h:38
TofFz_helix::r_etf
double r_etf
Definition Toffz_helix.h:42
TofFz_helix::Dz
double Dz
Definition Toffz_helix.h:38
TofFz_helix::Phi0
double Phi0
Definition Toffz_helix.h:38
TofFz_helix::Tofid
int Tofid
Definition Toffz_helix.h:26
TofFz_helix::W_tof
double W_tof
Definition Toffz_helix.h:30
TofFz_helix::Etfid
int Etfid
Definition Toffz_helix.h:27
TofFz_helix::Tanl
double Tanl
Definition Toffz_helix.h:38
TofFz_helix::R_tof
double R_tof
Definition Toffz_helix.h:28
TofFz_helix::Phi1
double Phi1
Definition Toffz_helix.h:40
TofFz_helix::Z_tof
double Z_tof
Definition Toffz_helix.h:34
TofFz_helix::Z_etf
double Z_etf
Definition Toffz_helix.h:35
TofFz_helix::Fi_tof
double Fi_tof
Definition Toffz_helix.h:29