CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DcxHel.cc
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // File and Version Information:
3 // $Id: DcxHel.cc,v 1.7 2012/01/31 16:41:41 wmtan Exp $
4 //
5 // Description:
6 // Class Implementation for |DcxHel|
7 //
8 // Environment:
9 // Software developed for the BaBar Detector at the SLAC B-Factory.
10 // Attempted port to CMSSW
11 //
12 // Author List:
13 // S. Wagner
14 //
15 //------------------------------------------------------------------------
16 //babar #include "BaBar/BaBar.hh"
17 //babar #include "BaBar/Constants.hh"
18 //babar #include "BbrGeom/BbrAngle.hh"
19 //babar #include <math.h>
20 //babar #include "DcxReco/DcxHel.hh"
21 //babar #include "DcxReco/DcxHit.hh"
22 #include <cmath>
23 #include "RecoTracker/RoadSearchHelixMaker/interface/DcxHel.hh"
24 #include "RecoTracker/RoadSearchHelixMaker/interface/DcxHit.hh"
25 
26 #include "FWCore/MessageLogger/interface/MessageLogger.h"
27 
28 using std::cout;
29 using std::endl;
30 
31 double DcxHel::half_pi=1.570796327;
32 double DcxHel::pi =3.141592654;
33 double DcxHel::twoPi =6.283185307;
34 //babar double DcxHel::ktopt =0.0045;
35 double DcxHel::ktopt =0.0120; //cms
36 
37 //constructors
38 
39 DcxHel::DcxHel( ){ }
40 
41 DcxHel::DcxHel
42 (double D0, double Phi0, double Omega, double Z0, double Tanl,
43  double T0, int Code, int Mode, double X, double Y){
44  omin=0.000005; omega=Omega; ominfl=1; if(fabs(omega)<omin){ominfl=0;}
45  phi0=Phi0; d0=D0;
46  if (phi0 > DcxHel::pi){phi0-=DcxHel::twoPi;}
47  if (phi0 < -DcxHel::pi){phi0+=DcxHel::twoPi;}
48  z0=Z0; tanl=Tanl;
49  xref=X; yref=Y; t0=T0;
50  cphi0=cos(phi0); sphi0=sin(phi0);
51  x0=X0(); y0=Y0(); xc=Xc(); yc=Yc();
52  code=Code;
53  decode(code,qd0,qphi0,qomega,qz0,qtanl,qt0,nfree);
54  mode=Mode; turnflag=0;
55 }//endof DcxHel
56 
57 DcxHel::DcxHel
58 (float D0, float Phi0, float Omega, float Z0, float Tanl,
59  float T0, int Code, int Mode, float X, float Y){
60  omin=0.000005; omega=Omega; ominfl=1; if(fabs(omega)<omin){ominfl=0;}
61  phi0=Phi0; d0=D0;
62  if (phi0 > DcxHel::pi){phi0-=DcxHel::twoPi;}
63  if (phi0 < -DcxHel::pi){phi0+=DcxHel::twoPi;}
64  z0=Z0; tanl=Tanl;
65  xref=X; yref=Y; t0=T0;
66  cphi0=cos(phi0); sphi0=sin(phi0);
67  x0=X0(); y0=Y0(); xc=Xc(); yc=Yc();
68  code=Code;
69  decode(code,qd0,qphi0,qomega,qz0,qtanl,qt0,nfree);
70  mode=Mode; turnflag=0;
71 }//endof DcxHel
72 
73 DcxHel::~DcxHel( ){ }
74 
75 //accessors
76 
77 double DcxHel::Xc()const{
78  if(ominfl){
79  return (X0()-sphi0/omega);
80  }else{
81  return 999999999.9;
82  }//(ominfl)
83 }//endof Xc
84 
85 double DcxHel::Yc()const{
86  if(ominfl){
87  return (Y0()+cphi0/omega);
88  }else{
89  return 999999999.9;
90  }//(ominfl)
91 }//endof Yc
92 
93 double DcxHel::X0()const{
94  return (xref-sphi0*d0);
95 }//endof X0
96 
97 double DcxHel::Y0()const{
98  return (yref+cphi0*d0);
99 }//endof Y0
100 
101 double DcxHel::Xh(double l)const{
102  if(ominfl){
103  double phit=phi0+omega*l;
104  return (xc+sin(phit)/omega);
105  }else{
106  return (x0+cphi0*l-0.5*l*l*omega*sphi0);
107  }//ominfl
108 }//endof Xh
109 
110 double DcxHel::Yh(double l)const{
111  if(ominfl){
112  double phit=phi0+omega*l;
113  return (yc-cos(phit)/omega);
114  }else{
115  return (y0+sphi0*l+0.5*l*l*omega*cphi0);
116  }//ominfl
117 }//endof Yh
118 
119 double DcxHel::Zh(double l)const{
120  return (z0+tanl*l);
121 }//endof Zh
122 
123 double DcxHel::Pt(double l)const{
124  if(ominfl){return DcxHel::ktopt/fabs(omega);}
125  else{return 1000000.0;}//ominfl
126 }//endof Px
127 
128 double DcxHel::Px(double l)const{
129  if(ominfl){double phit=phi0+omega*l; return DcxHel::ktopt*cos(phit)/fabs(omega);}
130  else{return 1000.0*cphi0;}//ominfl
131 }//endof Px
132 
133 double DcxHel::Py(double l)const{
134  if(ominfl){double phit=phi0+omega*l; return DcxHel::ktopt*sin(phit)/fabs(omega);}
135  else{return 1000.0*sphi0;}//ominfl
136 }//endof Py
137 
138 double DcxHel::Pz(double l)const{
139  if(ominfl){return DcxHel::ktopt*tanl/fabs(omega);}
140  else{return 1000.0*tanl;}//ominfl
141 }//endof Pz
142 
143 double DcxHel::Ptot(double l)const{
144  if(ominfl){return DcxHel::ktopt*sqrt(1.0+tanl*tanl)/fabs(omega);}
145  else{return 1000.0*sqrt(1.0+tanl*tanl);}//ominfl
146 }//endof Ptot
147 
148 double DcxHel::Lmax()const{
149  double lmax=250.0;
150  if(ominfl){
151  double rmax=1.0/fabs(omega);
152  double dmax=fabs(d0)+2.0*rmax;
153  if (dmax>80.0)lmax=DcxHel::pi*rmax;
154  }
155  return lmax;
156 }//endof Lmax
157 
158 //controls
159 //control fitting mode
160 void DcxHel::SetMode(int n){mode=n;}
161 void DcxHel::SetRef(double x, double y){xref=x; yref=y;}
162 //control free variables
163 void DcxHel::SetOmega(int Qomega){
164  nfree=nfree+deltaq(qomega,Qomega);
165  code=code+deltaq(qomega,Qomega)*100;
166  qomega=Qomega;
167 }
168 void DcxHel::SetPhi0(int Qphi0){
169  nfree=nfree+deltaq(qphi0,Qphi0);
170  code=code+deltaq(qphi0,Qphi0)*10;
171  qphi0=Qphi0;
172 }
173 void DcxHel::SetD0(int Qd0){
174  nfree=nfree+deltaq(qd0,Qd0);
175  code=code+deltaq(qd0,Qd0);
176  qd0=Qd0;
177 }
178 void DcxHel::SetTanl(int Qtanl){
179  nfree=nfree+deltaq(qtanl,Qtanl);
180  code=code+deltaq(qtanl,Qtanl)*10000;
181  qtanl=Qtanl;
182 }
183 void DcxHel::SetZ0(int Qz0){
184  nfree=nfree+deltaq(qz0,Qz0);
185  code=code+deltaq(qz0,Qz0)*1000;
186  qz0=Qz0;
187 }
188 void DcxHel::SetT0(int Qt0) {
189  nfree=nfree+deltaq(qt0, Qt0);
190  code=code+deltaq(qt0, Qt0)*100000;
191  qt0=Qt0;
192 }
193 
194 //operators
195 DcxHel& DcxHel::operator=(const DcxHel& rhs){
196  copy(rhs);
197  return *this;
198 }
199 
200 //decode free parameter code
201 void
202 DcxHel::decode(const int kode,int& i1,int& i2,
203  int& i3,int& i4,int& i5,int& i6,int& n)
204 {
205  int temp=kode;
206  temp=temp/1000000; temp=kode-1000000*temp;
207  i6=temp/100000; temp=temp-100000*i6;
208  i5=temp/10000; temp=temp-10000*i5;
209  i4=temp/1000; temp=temp-1000*i4;
210  i3=temp/100; temp=temp-100*i3;
211  i2=temp/10; i1=temp-10*i2;
212  n=0;
213  if(i6==1){n++;}else{i6=0;};
214  if(i5==1){n++;}else{i5=0;};
215  if(i4==1){n++;}else{i4=0;};
216  if(i3==1){n++;}else{i3=0;};
217  if(i2==1){n++;}else{i2=0;};
218  if(i1==1){n++;}else{i1=0;};
219 }//endof decode
220 
221 //copy |DcxHel| to |DcxHel|
222 void
223 DcxHel::copy(const DcxHel& rhs)
224 {
225  omega=rhs.Omega(); phi0=rhs.Phi0(); d0=rhs.D0(); t0=rhs.T0();
226  tanl=rhs.Tanl(); z0=rhs.Z0();
227  cphi0=rhs.CosPhi0(); sphi0=rhs.SinPhi0();
228  x0=rhs.X0(); y0=rhs.Y0(); xc=rhs.Xc(); yc=rhs.Yc();
229  xref=rhs.Xref(); yref=rhs.Yref();
230  qomega=rhs.Qomega(); qphi0=rhs.Qphi0(); qd0=rhs.Qd0(); qt0=rhs.Qt0();
231  qtanl=rhs.Qtanl(); qz0=rhs.Qz0();
232  mode=rhs.Mode(); nfree=rhs.Nfree();
233  code=rhs.Code(); ominfl=rhs.Ominfl(); omin=rhs.Omin();
234  turnflag=rhs.GetTurnFlag();
235 }//endof copy
236 
237 double
238 DcxHel::Doca( double wx, double wy, double wz,
239  double xi, double yi, double zi )
240 {
241  // describe wire
242 // edm::LogInfo("RoadSearch") << " In Doca, xi = " << xi << " yi = " << yi << " zi = " << zi ;
243  CLHEP::Hep3Vector ivec(xi,yi,zi);
244  wvec=CLHEP::Hep3Vector(wx,wy,wz);
245 // edm::LogInfo("RoadSearch") << " In Doca, wx = " << wx << " wy = " << wy << " wz = " << wz ;
246  // calculate len to doca
247  double xd=xi,yd=yi;
248 // edm::LogInfo("RoadSearch") << " In Doca, start xd = " << xd << " yd = " << yd ;
249  double lnew,t1,t2,dphi,dlen=1000.0; len=0.0; int itry=2;
250  // int segflg=0; if ((code==111)&&(z0==0.0)&&(tanl==0.0))segflg=1;
251  // int superseg=0; if ((code==11111)&&(xref!=0.0)&&(yref!=0.0))superseg=1;
252  double circut, circum=10000.;
253  if (ominfl){circum=DcxHel::twoPi/fabs(omega);} circut=0.50*circum;
254  while(itry){
255  if (ominfl){
256  t1=-xc+xd; t2=yc-yd; phi=atan2(t1,t2);
257  if (omega<0.0)phi+=DcxHel::pi; if (phi>DcxHel::pi)phi-=DcxHel::twoPi; dphi=phi-phi0;
258  if (omega < 0.0){
259  if (dphi > 0.0){dphi-=DcxHel::twoPi;}
260  if (dphi < -DcxHel::twoPi){dphi+=DcxHel::twoPi;}
261  }else{
262  if (dphi < 0.0){dphi+=DcxHel::twoPi;}
263  if (dphi > DcxHel::twoPi){dphi-=DcxHel::twoPi;}
264  }
265  lnew=dphi/omega;
266  // if ((lnew>circut)&&(segflg))lnew-=circum;
267  // if ((lnew>circut)&&(superseg))lnew-=circum;
268  if ((lnew>circut)&&(turnflag))lnew-=circum;
269  zh=Zh(lnew);
270  xd=xi+(zh-zi)*wx/wz;
271  yd=yi+(zh-zi)*wy/wz;
272  //double zd=zh;
273 // edm::LogInfo("RoadSearch") << " In Doca, xd = " << xd << " yd = " << yd << " zh = " << zh;
274 // edm::LogInfo("RoadSearch") << " lnew = " << lnew ;
275  dlen=fabs(lnew-len); len=lnew;
276  // if (segflg)break;
277  if (fabs(zh) > 250.0)break;
278  if ( (0.0==wx) && (0.0==wy) )break; if (dlen < 0.000001)break; itry--;
279  }else{len=(xi-xref)*cphi0+(yi-yref)*sphi0; zh=z0+tanl*len; phi=phi0; break;}
280  }
281  // CLHEP::Hep3Vector Dvec(xd,yd,zd);
282  xh=Xh(len); yh=Yh(len); CLHEP::Hep3Vector hvec(xh,yh,zh);
283 // edm::LogInfo("RoadSearch") << " In Doca, xh = " << xh << " yh = " << yh << " zh = " << zh ;
284  double lamb=atan(tanl); cosl=cos(lamb); sinl=sin(lamb);
285  tx=cosl*cos(phi); ty=cosl*sin(phi); tz=sinl;
286  tvec=CLHEP::Hep3Vector(tx,ty,tz);
287  CLHEP::Hep3Vector vvec=wvec.cross(tvec);
288  vhat=vvec.unit(); vx=vhat.x(); vy=vhat.y(); vz=vhat.z();
289 // edm::LogInfo("RoadSearch") << " In Doca, vx = " << vx << " vy = " << vy << " vz = " << vz ;
290  dvec=ivec-hvec; double doca=dvec*vhat;
291 // edm::LogInfo("RoadSearch") << " doca = " << doca ;
292  double f1=dvec*tvec; double f2=wvec*tvec; double f3=dvec*wvec;
293  f0=(f1-f2*f3)/(1.0-f2*f2);
294  if (doca>0.0){samb=-1;}else{samb=+1;}
295  double wirephi=atan2(yd,xd);
296  //babar eang=BbrAngle(phi-wirephi);
297  eang=normalize(phi-wirephi);
298  if (fabs(eang)<half_pi){wamb=samb;}else{wamb=-samb;}
299  if (fabs(zh) > 250.0)doca=1000.0;
300  return doca;
301 }//endof Doca
302 
303 std::vector<float>
304 DcxHel::derivatives(const DcxHit& hit)
305 {
306  double doca=Doca(hit);
307  std::vector<float> temp(nfree+1);
308  temp[0]=doca;
309  double fac=1.0;
310  if((mode==0)&&(doca<0.0)) {fac=-fac;}
311  if(mode==0) {temp[0]=fabs(temp[0]);}
312  int bump=0;
313  if(qd0){double dddd0=-vx*sphi0+vy*cphi0; bump++; temp[bump]=dddd0*fac;}
314  if(qphi0){
315  // double dddp0=-(yh-y0)*vx+(xh-x0)*vy;
316  double dddp0=-(yh-y0+f0*ty)*vx+(xh-x0+f0*tx)*vy;
317  dddp0=dddp0*(1.0+d0*omega);
318  bump++; temp[bump]=dddp0*fac;
319  }
320  if(qomega){double dddom;
321  if (ominfl){
322  dddom=((len*cos(phi)-xh+x0)*vx+(len*sin(phi)-yh+y0)*vy)/omega;
323  dddom+=f0*len*cosl*(-sin(phi)*vx+cos(phi)*vy);
324  }
325  else{dddom=0.5*len*len*(-vx*sphi0+vy*cphi0);}
326  bump++; temp[bump]=dddom*fac;
327  }
328  if(qz0){double dddz0=vz; bump++; temp[bump]=dddz0*fac;}
329  if(qtanl){double ddds=vz*len; bump++; temp[bump]=ddds*fac;}
330  if(qt0){bump++; temp[bump]=-hit.v();}
331  return temp;
332 }//endof derivatives
333 
334 void DcxHel::print()const {
335  edm::LogInfo("RoadSearch") << " "
336  << " d0 = " << d0
337  << " phi0 = " << phi0
338  << " omega = " << omega
339  << " z0 = " << z0
340  << " tanl = " << tanl
341  << " t0 = " << t0
342  << " code = " << code
343  << " mode = " << mode
344  << " ominfl = " << ominfl
345  << " nfree = " << nfree
346  << " x0 = " << x0
347  << " y0 = " << y0
348  << " xc = " << xc
349  << " yc = " << yc
350  << " xref = " << xref
351  << " yref = " << yref
352  << " " ;
353 }//endof print
354 
355 void DcxHel::flip(){
356  if (ominfl){
357  if ( (fabs(d0)+2.0/fabs(omega)) > 80.0)return;
358  double lturn=DcxHel::twoPi/fabs(omega); double zturn=Zh(lturn);
359 // edm::LogInfo("RoadSearch") << "z0 " << z0 << " zturn " << zturn ;
360  if (fabs(zturn) < fabs(z0)){
361  z0=zturn; tanl=-tanl; omega=-omega; d0=-d0;
362  phi0=phi0-DcxHel::pi; if (phi0<-DcxHel::pi){phi0+=DcxHel::twoPi;}
363  cphi0=cos(phi0); sphi0=sin(phi0); x0=X0(); y0=Y0();
364  }
365  }
366 }//endof flip
367 
368 double DcxHel::normalize(double angle) {
369  if (angle < - DcxHel::pi) {
370  angle += DcxHel::twoPi;
371  if (angle < - DcxHel::pi) angle = fmod(angle+ DcxHel::pi, DcxHel::twoPi) + DcxHel::pi;
372  }
373  else if (angle > DcxHel::pi) {
374  angle -= DcxHel::twoPi;
375  if (angle > DcxHel::pi) angle = fmod(angle+DcxHel::pi, DcxHel::twoPi) - DcxHel::pi;
376  }
377  return angle;
378 }
operator=
nocap nocap const skelname & operator=(const skelname &)
Z0
static const double Z0
Definition: CastorGeometryData.h:32
reco::print
std::string print(const Track &, edm::Verbosity=edm::Concise)
Track print utility.
Definition: print.cc:8
funct::D0
Divides< arg, void > D0
Definition: Factorize.h:143
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
X
#define X(str)
Definition: MuonsGrabber.cc:49
alignBH_cfg.mode
list mode
Definition: alignBH_cfg.py:7
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
edm::decode
bool decode(bool &, std::string const &)
Definition: types.cc:62
Geom::twoPi
double twoPi()
Definition: Pi.h:32
Y0
static const double Y0
Definition: CastorGeometryData.h:33
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
dmax
#define dmax(a, b)
Definition: mlp_lapack.h:164
pi
double pi
Definition: Vx3DHLTAnalyzer.h:63
x
Definition: DDAxes.h:10
X0
static const double X0
Definition: CastorGeometryData.h:31
angle
T angle(T x1, T y1, T z1, T x2, T y2, T z2)
Definition: angle.h:11
phi
Definition: DDAxes.h:10
Generated for CMSSW Reference Manual by doxygen 1.8.5