CGEM BOSS 6.6.5.g
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 _debug=1;
67 _pathl_cut= 500.0;
68 _ztof_cutm=-140.0;
69 _ztof_cutx= 140.0;
70 Tofid_mrpc=-1;
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 // fzisan track
96 /* FTFinder * ftrk = FTFinder:: GetPointer();
97 FTList<FTTrack *> & FsTrks = ftrk->tracks();
98 NTrk = FsTrks.length();
99
100 cout<< "Toffz_helix:: NTrk="<<NTrk<<endl;
101 if (NTrk<=0) return (-3);
102
103 // get helix params for Id_cdfztrk
104 FTTrack & trk = * FsTrks[Id_cdfztrk];
105 const Helix hel = * trk.helix();
106
107 const HepPoint3D pivot1(0.,0.,0.);
108 HepVector a(5,0);
109 a[0] = hel.a()[0];
110 a[1] = hel.a()[1];
111 a[2] = hel.a()[2];
112 a[3] = hel.a()[3];
113 a[4] = hel.a()[4];
114 */
115 const HepPoint3D pivot1(0.,0.,0.);
116 HepVector a(5,0);
117 a[0] = fzisan[0];
118 a[1] = fzisan[1];
119 a[2] = fzisan[2];
120 a[3] = fzisan[3];
121 a[4] = fzisan[4];
122
123
124 // Ill-fitted track in fzisan (dz=tanl=-9999.)
125 if (abs(a[3])>50.0 || abs(a[4])>500.0) return (-5);
126 // if (abs(a[3])>10.0 || abs(a[4])>500.0) return (-5);
127
128 // D e f i n e h e l i x
129 Helix helix1(pivot1,a);
130 Dr = helix1.a()[0];
131 Phi0 = helix1.a()[1];
132 Kappa = helix1.a()[2];
133 Dz = helix1.a()[3];
134 Tanl = helix1.a()[4];
135
136 /* double detaphi=asin(0.5*0.81*0.3*fabs(Kappa));
137 double phi=Phi0+detaphi+ piby1/2;
138 // int tofid = (int) ( phi/piby44 + 1.0 );
139 int tofid = (int) ( phi/piby44 );
140 */
141
142 // check cdfztrk hit on tof cyclinder
143 rho = helix1.radius();
144 // cout<<" Toffz_helix:: rho ="<<rho<<endl; // radius of the circle in cm
145 HepPoint3D xyz = helix1.center();
146 xc = xyz(0);
147 yc = xyz(1);
148 // cout<<"Toffz_helix::helix center:xc="<<xc<<",yc= "<<yc<<endl;
149
150 Hep3Vector vxyz = helix1.direction();
151 nx = vxyz(0); // direction of track at the vertex
152 ny = vxyz(1);
153 // cout<<"Toffz_helix::direction: nx= "<<nx<<", ny="<<ny<<endl;
154
155 if(fabs(rho)>R_tof/2){
156
157 // get hit on tof cylinder at radius R_tof;
158
159 if( xc==0.0 && yc==0.0 ) return(-3);
160 // coefficients of quadratic equation
161 ca = xc*xc + yc*yc ;
162 aa = 0.5 * ( ca - rho*rho + R_tof*R_tof );
163
164 if ( xc != 0.0 ) {
165 cb = aa * yc;
166 cc = aa*aa - R_tof*R_tof*xc*xc;
167 // determinant
168 detm = cb*cb - ca*cc;
169 if ( detm > 0.0 ) {
170 yy[0] = ( cb + sqrt(detm) )/ ca;
171 xx[0] = ( aa - yy[0]*yc )/xc;
172 yy[1] = ( cb - sqrt(detm) )/ ca;
173 xx[1] = ( aa - yy[1]*yc )/xc;
174 } else return(-1);
175 }
176 else{
177 cb = aa * xc;
178 cc = aa*aa - R_tof*R_tof*yc*yc;
179 // determinant
180 detm = cb*cb - ca*cc;
181 if ( detm > 0.0 ) {
182 xx[0] = ( cb + sqrt(detm) )/ca;
183 yy[0] = ( aa - xx[0]*xc )/yc;
184 xx[1] = ( cb - sqrt(detm) )/ca;
185 yy[1] = ( aa - xx[1]*xc )/yc;
186 } else return(-2);
187 }
188
189 // choose one of hits according to track direction at the vertex.
190
191 if( xx[0]*nx + yy[0]*ny > 0.0 ) { tofx = xx[0]; tofy = yy[0]; }
192 else { tofx = xx[1]; tofy = yy[1]; }
193
194 double fi = atan2(tofy ,tofx ); // atna2 from -pi to pi
195
196
197 if( fi < 0.0 ) fi=pi2+fi;
198 Fi_tof = fi;
199 // Tofid = (int) ( Fi_tof/piby44 + 1.0 );
200 Tofid = (int) ( Fi_tof/piby44 );
201
202 /* double fi2 = atan2(yy[1] , xx[1]); // atna2 from -pi to pi
203
204 if( fi2 < 0.0 ) fi2=pi2+fi2;
205 Fi_tof = fi2;
206 int Tofid2 = (int) ( Fi_tof/piby44 + 1.0 );
207
208 W_tof = ( Fi_tof - piby44*( (double) Tofid -1.0 ) )/ piby44; */
209
210 // get phi and z on the tof counter
211 const HepPoint3D pivot2( tofx, tofy, 0.0);
212 helix1.pivot(pivot2);
213
214 // corrected by Ozaki san to avoid negative path length on Aug.10,99
215 Phi1=helix1.a()[1];
216 Z_tof=helix1.a()[3];
217
218 double dphi = (-xc*(tofx-xc)-yc*(tofy-yc))/sqrt(xc*xc+yc*yc)/fabs(rho);
219 dphi = acos(dphi);
220 Pathl = fabs(rho*dphi);
221 double coscor=sqrt(1.0+Tanl*Tanl);
222 Pathl=Pathl*coscor;
223
224 // path length in tof scintillator
225 Hep3Vector vxyz1 = helix1.direction();
226 nx = vxyz1(0); // direction of track at the vertex
227 ny = vxyz1(1);
228 // cout<<" Toffz_helix::track direction: nx="<<nx<<",ny="<<ny<<endl;
229 double corxy=(nx*tofx+ny*tofy)/R_tof;
230 Path_tof=4.0*coscor/corxy;
231 // cout<<" Toffz_helix::Path_tof="<< Path_tof<<endl;
232 if(abs(Z_tof)<115) return (1);
233 // if(Z_tof>116.5){
234
235
236 //std::cout << " helix1.a()[3] " << helix1.a()[3] <<std::endl;
237
238 if(Z_tof>=115){
239 Z_tof=133;
240 Pathl=Z_tof/sin(atan(Tanl));
241 double phi0=-(Z_tof/Tanl)/rho;
242 double phi=-(Z_tof-Dz)/Tanl/rho;
243 double phi1 = (Z_tof*Kappa*0.003)/Tanl;
244 double phi2 = (Z_tof-Dz)*Kappa*0.003/Tanl;
245
246 double x_endtof=Dr*cos(Phi0)+(1/(Kappa*0.003))*(cos(Phi0)-cos(Phi0+phi));
247 double y_endtof=Dr*sin(Phi0)+(1/(Kappa*0.003))*(sin(Phi0)-sin(Phi0+phi));
248 r_endtof=sqrt(x_endtof * x_endtof + y_endtof * y_endtof);
249
250 Helix helix1(pivot1,a);
251
252 double x1_endtof=helix1.x(phi).x();
253 double y1_endtof=helix1.x(phi).y();
254 double z1_endtof=helix1.x(phi).z();
255
256 double fi_endtof = atan2(y_endtof,x_endtof ); // atna2 from -pi to pi
257 if (fi_endtof<0) fi_endtof=pi2+fi_endtof;
258
259 Tofid =(int)(fi_endtof/piby24);
260 if(Tofid>47)Tofid=Tofid-48;
261
262
263 double fi_endtof_mrpc = atan2(y_endtof,x_endtof)-3.141592654/2.; //For the MRPC detector the module number 1 is at x=0 y>0;
264 if(fi_endtof_mrpc<0) fi_endtof_mrpc += 2*3.141592654;
265
266 int mrpc_moduleId = ((int)(fi_endtof_mrpc/(3.141593/18)))+1; //We have 36 modules for the mrpc! My numbering starts with 1
267 if(mrpc_moduleId%2==0){mrpc_moduleId=mrpc_moduleId/2;}
268 else {mrpc_moduleId=(mrpc_moduleId+1)/2;}
269 if(helix1.a()[3]>0) {(mrpc_moduleId -= 19); mrpc_moduleId=mrpc_moduleId*(-1); }
270 Tofid_mrpc = mrpc_moduleId; //We only have an rough idea, what our module number could be! This has to be considerd later.
271 //std::cout << "Toffz helix Tofid_mrpc "<<Tofid_mrpc<<std::endl;
272
273
274 // }else if (Z_tof<-116.5){
275 return (0);
276
277 }else if (Z_tof<=-115){
278 Z_tof=-133;
279 Pathl=Z_tof/sin(atan(Tanl));
280 double phi0=-(Z_tof/Tanl)/rho;
281 double phi=-(Z_tof-Dz)/Tanl/rho;
282 double phi1 = (Z_tof*Kappa*0.003)/Tanl;
283 double phi2 = (Z_tof-Dz)*Kappa*0.003/Tanl;
284
285 double x_endtof=Dr*cos(Phi0)+(1/(Kappa*0.003))*(cos(Phi0)-cos(Phi0+phi));
286 double y_endtof=Dr*sin(Phi0)+(1/(Kappa*0.003))*(sin(Phi0)-sin(Phi0+phi));
287 r_endtof=sqrt(x_endtof * x_endtof + y_endtof * y_endtof);
288
289 Helix helix1(pivot1,a);
290
291 double x1_endtof=helix1.x(phi).x();
292 double y1_endtof=helix1.x(phi).y();
293 double z1_endtof=helix1.x(phi).z();
294
295 double fi_endtof = atan2(y_endtof,x_endtof ); // atna2 from -pi to pi
296 if (fi_endtof<0) fi_endtof=pi2+fi_endtof;
297
298 Tofid =(int)(fi_endtof/piby24);
299 if(Tofid>47) Tofid=Tofid-48;
300
301 double fi_endtof_mrpc = atan2(y_endtof,x_endtof)-3.141592654/2.; //For the MRPC detector the module number 1 is at x=0 y>0
302 if(fi_endtof_mrpc<0) fi_endtof_mrpc += 2*3.141592654;
303
304 int mrpc_moduleId = ((int)(fi_endtof_mrpc/(3.141593/18)))+1; //We have 36 modules for the mrpc! My numbering starts with 1
305 if(mrpc_moduleId%2==0){mrpc_moduleId=mrpc_moduleId/2;}
306 else {mrpc_moduleId=(mrpc_moduleId+1)/2;}
307 if(helix1.a()[3]>0) {(mrpc_moduleId -= 19); mrpc_moduleId=mrpc_moduleId*(-1); }
308 Tofid_mrpc = mrpc_moduleId; //We only have an rough idea, what our module number could be! This has to be considerd later.
309 //std::cout << "Toffz helix Tofid_mrpc "<< Tofid_mrpc<<std::endl;
310
311 return (2);
312
313
314 }
315
316 }
317 else {
318
319 if(Tanl>0){
320 Z_tof=133;
321
322 Pathl=Z_tof/sin(atan(Tanl));
323 double phi0=-(Z_tof/Tanl)/rho;
324 double phi=-(Z_tof-Dz)/Tanl/rho;
325 double phi1 = (Z_tof*Kappa*0.003)/Tanl;
326 double phi2 = (Z_tof-Dz)*Kappa*0.003/Tanl;
327
328 double x_endtof=Dr*cos(Phi0)+(1/(Kappa*0.003))*(cos(Phi0)-cos(Phi0+phi));
329 double y_endtof=Dr*sin(Phi0)+(1/(Kappa*0.003))*(sin(Phi0)-sin(Phi0+phi));
330 r_endtof=sqrt(x_endtof * x_endtof + y_endtof * y_endtof);
331 // double x2_endtof=Dr*cos(Phi0)+(1/(Kappa*0.003))*(cos(Phi0)-cos(Phi0+phi2));
332 //double y2_endtof=Dr*sin(Phi0)+(1/(Kappa*0.003))*(sin(Phi0)-sin(Phi0+phi2));
333 //double r2_endtof=sqrt(x2_endtof * x2_endtof + y2_endtof * y2_endtof);
334
335 double fi_endtof = atan2(y_endtof,x_endtof ); // atna2 from -pi to pi
336 if (fi_endtof<0) fi_endtof=pi2+fi_endtof;
337
338 Tofid =(int)(fi_endtof/piby24);
339 if(Tofid>47)Tofid=Tofid-48;
340
341
342 double fi_endtof_mrpc = atan2(y_endtof,x_endtof)-3.141592654/2.; //For the MRPC detector the module number 1 is at x=0 y>0
343 if(fi_endtof_mrpc<0) fi_endtof_mrpc += 2*3.141592654;
344
345 int mrpc_moduleId = ((int)(fi_endtof_mrpc/(3.141593/18)))+1; //We have 36 modules for the mrpc! My numbering starts with 1
346 if(mrpc_moduleId%2==0){mrpc_moduleId=mrpc_moduleId/2;}
347 else {mrpc_moduleId=(mrpc_moduleId+1)/2;}
348 if(helix1.a()[3]>0) {(mrpc_moduleId -= 19); mrpc_moduleId=mrpc_moduleId*(-1); }
349 Tofid_mrpc = mrpc_moduleId; //We only have an rough idea, what our module number could be! This has to be considerd later.
350 //std::cout << "Toffz helix Tofid_mrpc "<< Tofid_mrpc<<std::endl;
351
352
353
354
355 return (0);
356
357 }
358 else{
359 Z_tof=-133;
360 Pathl=Z_tof/sin(atan(Tanl));
361 double phi0=-(Z_tof/Tanl)/rho;
362 double phi=-(Z_tof-Dz)/Tanl/rho;
363 double phi1 = (Z_tof*Kappa*0.003)/Tanl;
364 double phi2 = (Z_tof-Dz)*Kappa*0.003/Tanl;
365
366
367 double x_endtof=Dr*cos(Phi0)+(1/(Kappa*0.003))*(cos(Phi0)-cos(Phi0+phi));
368 double y_endtof=Dr*sin(Phi0)+(1/(Kappa*0.003))*(sin(Phi0)-sin(Phi0+phi));
369 r_endtof=sqrt(x_endtof * x_endtof + y_endtof * y_endtof);
370 //double x2_endtof=Dr*cos(Phi0)+(1/(Kappa*0.003))*(cos(Phi0)-cos(Phi0+phi2));
371 //double y2_endtof=Dr*sin(Phi0)+(1/(Kappa*0.003))*(sin(Phi0)-sin(Phi0+phi2));
372 //double r2_endtof=sqrt(x2_endtof * x2_endtof + y2_endtof * y2_endtof);
373
374 double fi_endtof = atan2(y_endtof,x_endtof ); // atna2 from -pi to pi
375 if (fi_endtof<0) fi_endtof=pi2+fi_endtof;
376
377 Tofid =(int)(fi_endtof/piby24);
378 if(Tofid>47)Tofid=Tofid-48;
379
380
381 double fi_endtof_mrpc = atan2(y_endtof,x_endtof)-3.141592654/2.; //For the MRPC detector the module number 1 is at x=0 y>0
382 if(fi_endtof_mrpc<0) fi_endtof_mrpc += 2*3.141592654;
383
384 int mrpc_moduleId = ((int)(fi_endtof_mrpc/(3.141593/18)))+1; //We have 36 modules for the mrpc! My numbering starts with 1
385 if(mrpc_moduleId%2==0){mrpc_moduleId=mrpc_moduleId/2;}
386 else {mrpc_moduleId=(mrpc_moduleId+1)/2;}
387 if(helix1.a()[3]>0) {(mrpc_moduleId -= 19); mrpc_moduleId=mrpc_moduleId*(-1); }
388 Tofid_mrpc = mrpc_moduleId; //We only have an rough idea, what our module number could be! This has to be considerd later.
389 //std::cout << "Toffz helix Tofid_mrpc "<< Tofid_mrpc<<std::endl;
390
391
392
393
394
395 return (2);
396 }
397 }
398
399 if (abs(Pathl) > _pathl_cut ||
400 Z_tof < _ztof_cutm || Z_tof > _ztof_cutx) {
401 // Bad path length or Z_tof
402 if (_debug) {
403 printf("\n TofFz_helix> Trk=%3d params(dr,phi,kappa,dz,tanl)="
404 "(%5.1f,%6.3f,%6.4f,%6.1f,%6.3f) R_tof %5.1f\n",
405 Id_cdfztrk,Dr,Phi0,Kappa,Dz,Tanl,R_tof);
406
407 printf(" TofFz_helix> rho=%8.1f, (xc,yc)=(%8.1f,%8.1f)"
408 " (nx,ny)=(%5.2f,%5.2f)\n",rho,xc,yc,nx,ny);
409
410 printf(" TofFz_helix> tof (x,y)=(%5.1f,%5.1f) and (%5.1f,%5.1f)\n",
411 xx[0],yy[0],xx[1],yy[1]);
412
413 printf(" TofFz_helix> tofid=%3d, fitof=%6.3f, w=%5.3f"
414 " (x,y,z)=(%5.1f,%5.1f,%5.1f) pathl=%5.1f cm path=%5.1f cm\n",
415 Tofid,Fi_tof,W_tof,tofx,tofy,Z_tof,Pathl,Path_tof);
416 }
417 return (-7);
418 }
419
420
421 // return(1);
422
423}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
phi2
Double_t phi2
Definition: DataBase/tau_mode.c:7
phi1
Double_t phi1
Definition: DataBase/tau_mode.c:7
abs
double abs(const EvtComplex &c)
Definition: EvtComplex.hh:212
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-12/TrackUtil/Helix.h:212
Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/src/Helix.cxx:152
Helix::radius
double radius(void) const
returns radious of helix.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:224
Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:272
Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:218
Helix::direction
Hep3Vector direction(double dPhi=0.) const
returns direction vector after rotating angle dPhi in phi direction.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:230
TofFz_helix::TofFz_Get
int TofFz_Get(double, int, double[])
Definition: Toffz_helix.cxx:76
TofFz_helix::Kappa
double Kappa
Definition: Toffz_helix.h:37
TofFz_helix::Path_tof
double Path_tof
Definition: Toffz_helix.h:33
TofFz_helix::Tofid_mrpc
int Tofid_mrpc
Definition: Toffz_helix.h:28
TofFz_helix::r_endtof
double r_endtof
Definition: Toffz_helix.h:40
TofFz_helix::Pathl
double Pathl
Definition: Toffz_helix.h:32
TofFz_helix::TofFz_helix
TofFz_helix(void)
Definition: Toffz_helix.cxx:61
TofFz_helix::Dr
double Dr
Definition: Toffz_helix.h:37
TofFz_helix::Dz
double Dz
Definition: Toffz_helix.h:37
TofFz_helix::Phi0
double Phi0
Definition: Toffz_helix.h:37
TofFz_helix::Tofid
int Tofid
Definition: Toffz_helix.h:27
TofFz_helix::W_tof
double W_tof
Definition: Toffz_helix.h:31
TofFz_helix::Tanl
double Tanl
Definition: Toffz_helix.h:37
TofFz_helix::R_tof
double R_tof
Definition: Toffz_helix.h:29
TofFz_helix::Phi1
double Phi1
Definition: Toffz_helix.h:39
TofFz_helix::Z_tof
double Z_tof
Definition: Toffz_helix.h:34
TofFz_helix::Fi_tof
double Fi_tof
Definition: Toffz_helix.h:30