00001 #include "MdcUtilitySvc/MdcUtilitySvc.h"
00002 #include "GaudiKernel/MsgStream.h"
00003 #include "GaudiKernel/SmartDataPtr.h"
00004 #include "GaudiKernel/Bootstrap.h"
00005 #include "GaudiKernel/ISvcLocator.h"
00006 #include "EventModel/EventHeader.h"
00007 #include "MdcGeom/Constants.h"
00008 #include "MdcGeom/BesAngle.h"
00009 #include "TrkBase/HelixTraj.h"
00010 #include "TrkBase/TrkPoca.h"
00011 #include "MdcGeom/MdcLayer.h"
00012 #include "GaudiKernel/IDataProviderSvc.h"
00013 #include <cmath>
00014 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00015
00016 typedef HepGeom::Point3D<double> HepPoint3D;
00017 #endif
00018
00019 using namespace std;
00020 using namespace Event;
00021
00022 MdcUtilitySvc::MdcUtilitySvc( const string& name, ISvcLocator* svcloc) :
00023 Service (name, svcloc) {
00024 declareProperty("debug", m_debug = 0);
00025 }
00026
00027 MdcUtilitySvc::~MdcUtilitySvc(){
00028 }
00029
00030 StatusCode MdcUtilitySvc::initialize(){
00031 MsgStream log(messageService(), name());
00032 log << MSG::INFO << "MdcUtilitySvc::initialize()" << endreq;
00033
00034 StatusCode sc = Service::initialize();
00035 if( sc.isFailure() ) {
00036 log << MSG::ERROR << "Service::initialize() failure" << endreq;
00037 return sc;
00038 }
00039
00040
00041 sc = service("MagneticFieldSvc", m_pIMF);
00042 if(! m_pIMF){
00043 log << MSG::FATAL <<" ERROR Unable to open Magnetic field service "<< endreq;
00044 return StatusCode::FAILURE;
00045 }
00046
00047
00048 sc = service("MdcGeomSvc", m_mdcGeomSvc);
00049 if(! m_mdcGeomSvc){
00050 log << MSG::FATAL <<" Could not load MdcGeomSvc! "<< endreq;
00051 return StatusCode::FAILURE;
00052 }
00053
00054 return StatusCode::SUCCESS;
00055 }
00056
00057 StatusCode MdcUtilitySvc::finalize(){
00058 MsgStream log(messageService(), name());
00059 log << MSG::INFO << "MdcUtilitySvc::finalize()" << endreq;
00060
00061 return StatusCode::SUCCESS;
00062 }
00063
00064 StatusCode MdcUtilitySvc::queryInterface(const InterfaceID& riid, void** ppvInterface){
00065
00066 if( IID_IMdcUtilitySvc.versionMatch(riid) ){
00067 *ppvInterface = static_cast<IMdcUtilitySvc*> (this);
00068 } else{
00069 return Service::queryInterface(riid, ppvInterface);
00070 }
00071
00072 addRef();
00073 return StatusCode::SUCCESS;
00074 }
00075
00076
00077 int MdcUtilitySvc::nLayerTrackPassed(const double helixBes[5]) const{
00078
00079 HepVector helixBesParam(5,0);
00080 for(int i=0; i<5; ++i) helixBesParam[i] = helixBes[i];
00081
00082 return nLayerTrackPassed(helixBesParam);
00083 }
00084
00085
00086 int MdcUtilitySvc::nLayerTrackPassed(const HepVector helixBes) const{
00087 HepVector helix= besPar2PatPar(helixBes);
00088
00089 int nLayer = 0;
00090
00091 for(unsigned iLayer=0; iLayer<43; iLayer++){
00092
00093
00094 double rMidLayer = m_mdcGeomSvc->Layer(iLayer)->Radius();
00095 double flightLength = rMidLayer;
00096
00097 HepPoint3D pivot(0.,0.,0.);
00098 double dz = helix[3];
00099 double c = CLHEP::c_light * 100.;
00100 double alpha = 1/(c * Bz());
00101 double kappa = helix[2];
00102 double rc = (-1.)*alpha/kappa;
00103
00104 double tanl = helix[4];
00105 double phi0 = helix[1];
00106 double phi = flightLength/rc + phi0;
00107 double z = pivot.z() + dz - (alpha/kappa) * tanl * phi;
00108
00109 double layerHalfLength = m_mdcGeomSvc->Layer(iLayer)->Length()/2.;
00110
00111
00112
00113 if (fabs(z) < fabs(layerHalfLength)) ++nLayer;
00114 }
00115
00116 return nLayer;
00117 }
00118
00119
00120 HepVector MdcUtilitySvc::patPar2BesPar(const HepVector& helixPar) const{
00121 HepVector helix(5,0);
00122 double d0 = -helixPar[0];
00123 double phi0 = helixPar[1]+ CLHEP::halfpi;
00124 double omega = Bz()*helixPar[2]/-333.567;
00125 double z0 = helixPar[3];
00126 double tanl = helixPar[4];
00127 helix[0] = d0;
00128 helix[1] = phi0;
00129 helix[2] = omega;
00130 helix[3] = z0;
00131 helix[4] = tanl;
00132
00133 return helix;
00134 }
00135
00136
00137 HepSymMatrix MdcUtilitySvc::patErr2BesErr(const HepSymMatrix& err) const{
00138
00139
00140
00141
00142 HepSymMatrix mS(err.num_row(),0);
00143 mS[0][0]=-1.;
00144 mS[1][1]=1.;
00145
00146 mS[2][2]=Bz()/-333.567;
00147 mS[3][3]=1.;
00148 mS[4][4]=1.;
00149 HepSymMatrix mVy= err.similarity(mS);
00150
00151 return mVy;
00152 }
00153
00154
00155
00156 HepPoint3D MdcUtilitySvc::pointOnHelix(const HepVector helixBes, int layer, int innerOrOuter) const{
00157 HepVector helixPat= besPar2PatPar(helixBes);
00158 return pointOnHelixPatPar(helixPat, layer, innerOrOuter);
00159 }
00160
00161
00162
00163 HepPoint3D MdcUtilitySvc::pointOnHelixPatPar(const HepVector helixPat, int layer, int innerOrOuter) const{
00164
00165 double rInner,rOuter;
00166
00167 double rCize1 =0.1* m_mdcGeomSvc->Layer(layer)->RCSiz1();
00168 double rCize2 =0.1* m_mdcGeomSvc->Layer(layer)->RCSiz2();
00169 double rLay =0.1* m_mdcGeomSvc->Layer(layer)->Radius();
00170
00171
00172 rInner = rLay - rCize1 ;
00173 rOuter = rLay + rCize2 ;
00174
00175
00176 HepPoint3D piv(0.,0.,0.);
00177 double r;
00178 if (innerOrOuter) r = rInner;
00179 else r = rOuter;
00180
00181 double d0 = helixPat[0];
00182 double phi0 = helixPat[1];
00183 double omega = helixPat[2];
00184 double z0 = helixPat[3];
00185 double tanl = helixPat[4];
00186
00187 double rc;
00188 if (abs(omega) <Constants::epsilon) rc = 9999.;
00189 else rc = 1./omega;
00190 double xc = piv.x() + ( d0 + rc) * cos(phi0);
00191 double yc = piv.y() + ( d0 + rc) * sin(phi0);
00192 HepPoint3D helixCenter(xc,yc,0.0);
00193 rc = sqrt(xc*xc + yc*yc);
00194 double a,b,c;
00195 a = r;
00196 b = d0 + rc;
00197 c = rc;
00198 double dphi = acos((a*a-b*b-c*c)/(-2.*b*c));
00199 double fltlen = dphi * rc;
00200 double phi = phi0 * omega * fltlen;
00201 double x = piv.x()+d0*sin(phi) - (rc+d0)*sin(phi0);
00202 double y = piv.y()+d0*cos(phi) + (rc+d0)*cos(phi0);
00203 double z = piv.z()+ z0 + fltlen*tanl;
00204
00205 if(m_debug) {
00206 std::cout<<" abc "<<a<<" "<<b<<" "<<c<<" omega "<<omega<<" r "<<r<<" rc "<<rc<<" dphi "<<dphi<<" piv "<<piv.z()
00207 <<" z0 "<<z0<<" fltlen "<<fltlen<<" tanl "<<tanl<<std::endl;
00208 std::cout<<" pointOnHelixPatPar in Hel "<<helixPat<<std::endl;
00209 cout<<"HepPoint3D(x, y, z) = "<<HepPoint3D(x, y, z)<<endl;
00210 }
00211 return HepPoint3D(x, y, z);
00212 }
00213
00214
00215
00216 double MdcUtilitySvc::doca(int layer, int cell, const HepVector helixBes, const HepSymMatrix errMatBes, bool passCellRequired, bool doSag) const{
00217 HepVector helixPat = besPar2PatPar(helixBes);
00218 HepSymMatrix errMatPat= besErr2PatErr(errMatBes);
00219
00220 return docaPatPar(layer, cell, helixPat, errMatPat, passCellRequired, doSag);
00221 }
00222
00223
00224
00225 double MdcUtilitySvc::doca(int layer, int cell,HepPoint3D eastP,HepPoint3D westP, const HepVector helixBes, const HepSymMatrix errMatBes, bool passCellRequired, bool doSag) const{
00226 HepVector helixPat = besPar2PatPar(helixBes);
00227 HepSymMatrix errMatPat= besErr2PatErr(errMatBes);
00228
00229 return docaPatPar(layer, cell, eastP, westP, helixPat, errMatPat, passCellRequired, doSag);
00230 }
00231
00232
00233
00234 double MdcUtilitySvc::doca(int layer, int cell, const MdcSWire* sWire, const HepVector helixBes,const HepSymMatrix errMatBes, bool passCellRequired) const{
00235 HepVector helixPat = besPar2PatPar(helixBes);
00236 HepSymMatrix errMatPat= besErr2PatErr(errMatBes);
00237
00238 return docaPatPar(layer, cell, sWire, helixPat, errMatPat, passCellRequired);
00239 }
00240
00241
00242
00243 double MdcUtilitySvc::docaPatPar(int layer, int cell, const HepVector helixPat, const HepSymMatrix errMat, bool passCellRequired, bool doSag) const{
00244
00245 const MdcGeoWire *geowir = m_mdcGeomSvc->Wire(layer,cell);
00246 int id = geowir->Id();
00247 double sag = 0.;
00248 if(doSag) sag = geowir->Sag()/10.;
00249 HepPoint3D eastP = geowir->Backward()/10.0;
00250 HepPoint3D westP = geowir->Forward() /10.0;
00251 const MdcSWire* sWire = new MdcSWire(eastP, westP, sag, id , cell);
00252
00253 double doca = docaPatPar(layer, cell, sWire, helixPat, errMat, passCellRequired);
00254
00255 delete sWire;
00256 return doca;
00257 }
00258
00259
00260
00261 double MdcUtilitySvc::docaPatPar(int layer, int cell,HepPoint3D eastP,HepPoint3D westP, const HepVector helixPat,const HepSymMatrix errMat, bool passCellRequired, bool doSag) const{
00262
00263 const MdcGeoWire *geowir = m_mdcGeomSvc->Wire(layer,cell);
00264 int id = geowir->Id();
00265 double sag = 0.;
00266 if(doSag) sag = geowir->Sag()/10.;
00267 const MdcSWire* sWire = new MdcSWire(eastP, westP, sag, id , cell);
00268
00269 double doca = docaPatPar(layer, cell, sWire, helixPat, errMat, passCellRequired);
00270
00271 delete sWire;
00272 return doca;
00273 }
00274
00275
00276
00277 double MdcUtilitySvc::docaPatPar(int layer, int cell, const MdcSWire* sWire, const HepVector helixPat,const HepSymMatrix errMat, bool passCellRequired) const{
00278
00279 if(m_debug) std::cout<<" helixPat "<<helixPat<<std::endl;
00280 if(m_debug) std::cout<<" layer "<<layer<<" cell "<<cell<<std::endl;
00281
00282 int cellId_in = -1;
00283 int cellId_out = -1;
00284
00285
00286 cellTrackPassedByPhiPatPar(helixPat,layer,cellId_in,cellId_out);
00287
00288 if(m_debug) {
00289 std::cout<<" cellId in "<<cellId_in<<" out "<<cellId_out <<std::endl;
00290 cout << "cell = " << cell << ", cellId_in = " << cellId_in << ", cellId_out = " << cellId_out << endl;
00291 }
00292 if (passCellRequired &&(cell < cellId_in && cell > cellId_out)) return -999.;
00293
00294
00295 HelixTraj* helixTraj = new HelixTraj(helixPat,errMat);
00296 const Trajectory* trajHelix = dynamic_cast<const Trajectory*> (helixTraj);
00297
00298
00299 int innerOrOuter = 1;
00300 Hep3Vector cell_IO = pointOnHelixPatPar(helixPat,layer,innerOrOuter);
00301 double fltTrack = 0.1 * m_mdcGeomSvc -> Layer(layer)->Radius();
00302 if(m_debug) {
00303 std::cout<<" cell_IO "<<cell_IO<<std::endl;
00304 std::cout<<" fltTrack "<<fltTrack<<std::endl;
00305 }
00306
00307
00308 const MdcSagTraj* m_wireTraj = sWire->getTraj();
00309 const Trajectory* trajWire = dynamic_cast<const Trajectory*> (m_wireTraj);
00310 const HepPoint3D* start_In = sWire->getEastPoint();
00311
00312 if(m_debug){
00313 std::cout<<" sag "<<m_wireTraj->sag()<<std::endl;
00314 std::cout<< " east -------- "<< start_In->x()<<","<<start_In->y()<<","<<start_In->z()<<std::endl;
00315 }
00316
00317
00318
00319 double doca = -999.;
00320 TrkPoca* trkPoca;
00321 double zWire = cell_IO.z();
00322
00323 HepPoint3D pos_in(sWire->xWireDC(zWire),sWire->yWireDC(zWire),zWire) ;
00324 if(m_debug) std::cout<<" zWire "<<zWire<<" zEndDC "<<sWire->zEndDC()<<std::endl;
00325
00326
00327 double fltWire = sqrt( (pos_in.x()-start_In->x())*(pos_in.x()-start_In->x()) +
00328 (pos_in.y()-start_In->y())*(pos_in.y()-start_In->y()) +
00329 (pos_in.z()-start_In->z())*(pos_in.z()-start_In->z()) );
00330 trkPoca = new TrkPoca(*trajHelix, fltTrack, *trajWire, fltWire);
00331 doca = trkPoca->doca();
00332
00333 double hitLen = trkPoca->flt2();
00334 double startLen = trajWire->lowRange() - 5.;
00335 double endLen = trajWire->hiRange() + 5.;
00336 if(hitLen < startLen || hitLen > endLen) {
00337 if(m_debug) std::cout<<"WARNING MdcUtilitySvc::docaPatPar() isBeyondEndflange! hitLen="<<hitLen <<" startLen="<<startLen<<" endLen "<<endLen<<std::endl;
00338 doca = -998.;
00339 }
00340
00341
00342 if(m_debug) std::cout<<" doca "<<doca<<std::endl;
00343 delete helixTraj;
00344 delete trkPoca;
00345
00346 return doca;
00347
00348 }
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378 bool MdcUtilitySvc::cellTrackPassedByPhi(const HepVector helixBes,int layer,int& cellId_in,int& cellId_out) const{
00379 HepVector helixPat = besPar2PatPar(helixBes);
00380 return cellTrackPassedByPhiPatPar(helixPat, layer, cellId_in, cellId_out);
00381 }
00382
00383
00384
00385
00386
00387 bool MdcUtilitySvc::cellTrackPassedByPhiPatPar(const HepVector helixPat,int layer,int& cellId_in,int& cellId_out) const{
00388 int nCell = m_mdcGeomSvc->Layer(layer)->NCell();
00389
00390 double dPhiz = (m_mdcGeomSvc->Wire(layer,0)->Forward().phi() - m_mdcGeomSvc->Wire(layer,0)->Backward().phi())*0.5;
00391 double rEnd = m_mdcGeomSvc->Wire(layer,0)->Backward().rho()/10.;
00392 double rMid = rEnd * cos(dPhiz);
00393
00394
00395 double fltLenIn = rMid;
00396 double phiIn = helixPat[1] + helixPat[2]*fltLenIn;
00397
00398 double phiEPOffset = m_mdcGeomSvc->Wire(layer,0)->Backward().phi();
00399 BesAngle tmp(phiIn - phiEPOffset);
00400 if(m_debug) std::cout<<"cellTrackPassed nCell "<<nCell<<" layer "<<layer<<" fltLenIn "<<fltLenIn<<" phiEPOffset "<<phiEPOffset<<std::endl;
00401
00402 int wlow = (int)floor(nCell * tmp.rad() / twopi );
00403 int wbig = (int)ceil(nCell * tmp.rad() / twopi );
00404 if (tmp == 0 ){ wlow = -1; wbig = 1; }
00405
00406 if ((wlow%nCell)< 0) wlow += nCell;
00407 cellId_in = wlow;
00408
00409 if ((wbig%nCell)< 0) wbig += nCell;
00410 cellId_out = wbig;
00411
00412 bool passedOneCell = (cellId_in == cellId_out);
00413
00414 return passedOneCell;
00415 }
00416
00417
00418 bool MdcUtilitySvc::cellTrackPassedPatPar(HepVector helixPat, int layer,int& cellId_in,int& cellId_out) const{
00419 HepVector helixBes = patPar2BesPar(helixPat);
00420 return cellTrackPassed(helixBes, layer, cellId_in, cellId_out);
00421 }
00422
00423
00424
00425 bool MdcUtilitySvc::cellTrackPassed(HepVector helixBes, int lay,int& cellId_in,int& cellId_out) const{
00426 int charge,type,nCell;
00427 double dr0,phi0,kappa,dz0,tanl;
00428 double ALPHA_loc,rho,phi,cosphi0,sinphi0,x_hc,y_hc,z_hc;
00429 double dphi0,IO_phi,C_alpha,xx,yy;
00430 double inlow,inup,outlow,outup,phi_in,phi_out,phi_bin;
00431 double rCize1,rCize2,rLay,length,phioffset,slant,shift;
00432 double m_crio[2], phi_io[2], stphi[2],phioff[2],dphi[2];
00433
00434 dr0 = helixBes[0];
00435 phi0 = helixBes[1];
00436 kappa = helixBes[2];
00437 dz0 = helixBes[3];
00438 tanl = helixBes[4];
00439
00440 ALPHA_loc = 1000/(2.99792458*Bz());
00441 charge = ( kappa >=0 ) ? 1 : -1 ;
00442 rho = ALPHA_loc/kappa ;
00443 double pi = Constants::pi;
00444 phi = fmod(phi0 + 4*pi , 2*pi);
00445 cosphi0 = cos(phi);
00446 sinphi0 = (1.0 - cosphi0 ) * (1.0 + cosphi0 );
00447 sinphi0 = sqrt(( sinphi0 > 0.) ? sinphi0 : 0.);
00448 if( phi > pi ) sinphi0 = -sinphi0 ;
00449
00450 x_hc = 0. + ( dr0 + rho ) * cosphi0;
00451 y_hc = 0. + ( dr0 + rho ) * sinphi0;
00452 z_hc = 0. + dz0;
00453
00454
00455 HepPoint3D piv(0.,0.,0.);
00456 HepPoint3D hcenter(x_hc,y_hc,0.0);
00457
00458 double m_c_perp(hcenter.perp());
00459 Hep3Vector m_c_unit((HepPoint3D)hcenter.unit());
00460 HepPoint3D IO = (-1) * charge * m_c_unit;
00461
00462 dphi0 = fmod(IO.phi()+4*pi, 2*pi) - phi;
00463 IO_phi = fmod(IO.phi()+4*pi, 2*pi);
00464
00465 if(dphi0 > pi) dphi0 -= 2*pi;
00466 else if(dphi0 < -pi) dphi0 += 2*pi;
00467
00468 phi_io[0] = -(1+charge)*pi/2 - charge*dphi0;
00469 phi_io[1] = phi_io[0]+1.5*pi;
00470
00471
00472 bool outFlag = false;
00473
00474 rCize1 = 0.1 * m_mdcGeomSvc->Layer(lay)->RCSiz1();
00475 rCize2 = 0.1 * m_mdcGeomSvc->Layer(lay)->RCSiz2();
00476 rLay = 0.1 * m_mdcGeomSvc->Layer(lay)->Radius();
00477 length = 0.1 * m_mdcGeomSvc->Layer(lay)->Length();
00478
00479
00480 nCell = m_mdcGeomSvc->Layer(lay)->NCell();
00481 phioffset = m_mdcGeomSvc->Layer(lay)->Offset();
00482 slant = m_mdcGeomSvc->Layer(lay)->Slant();
00483 shift = m_mdcGeomSvc->Layer(lay)->Shift();
00484 type = m_mdcGeomSvc->Layer(lay)->Sup()->Type();
00485
00486 m_crio[0] = rLay - rCize1 ;
00487 m_crio[1] = rLay + rCize2 ;
00488
00489 int sign = -1;
00490
00491 Hep3Vector iocand[2];
00492 Hep3Vector cell_IO[2];
00493
00494 for(int ii =0; ii<2; ii++){
00495
00496 double cos_alpha = (m_c_perp*m_c_perp + m_crio[ii]*m_crio[ii] - rho*rho)
00497 / ( 2 * m_c_perp * m_crio[ii] );
00498 if(fabs(cos_alpha)>1&&(ii==0)){
00499 cos_alpha = -1;
00500 outFlag=true;
00501 }
00502 if(fabs(cos_alpha)>1&&(ii==1)){
00503 cos_alpha = (m_c_perp*m_c_perp + m_crio[0]*m_crio[0] - rho*rho)
00504 / ( 2 * m_c_perp * m_crio[0] );
00505 C_alpha = 2*pi - acos( cos_alpha);
00506 }else {
00507 C_alpha = acos( cos_alpha );
00508 }
00509
00510 iocand[ii] = m_c_unit;
00511 iocand[ii].rotateZ( charge*sign*C_alpha );
00512 iocand[ii] *= m_crio[ii];
00513
00514 xx = iocand[ii].x() - x_hc ;
00515 yy = iocand[ii].y() - y_hc ;
00516
00517 dphi[ii] = atan2(yy,xx) - phi0 - 0.5*pi*(1-charge);
00518 dphi[ii] = fmod( dphi[ii] + 8.0*pi,2*pi);
00519
00520
00521 if( dphi[ii] < phi_io[0] ) {
00522 dphi[ii] += 2*pi;
00523 }else if( dphi[ii] > phi_io[1] ){
00524 dphi[ii] -= 2*pi;
00525 }
00526
00527 cell_IO[ii] = Hel(piv, dr0, phi, ALPHA_loc, kappa,dz0, dphi[ii], tanl);
00528
00529
00530 if( (cell_IO[ii].x()==0 ) && (cell_IO[ii].y()==0) && (cell_IO[ii].z()==0)) continue;
00531 }
00532
00533
00534 cellId_in = cellId_out = -1 ;
00535 phi_in = cell_IO[0].phi();
00536 phi_in = fmod ( phi_in + 4 * pi, 2 * pi );
00537 phi_out = cell_IO[1].phi();
00538 phi_out = fmod ( phi_out + 4 * pi, 2 * pi );
00539 phi_bin = 2.0 * pi / nCell ;
00540
00541
00542 stphi[0] = shift * phi_bin * (0.5 - cell_IO[0].z()/length);
00543 stphi[1] = shift * phi_bin * (0.5 - cell_IO[1].z()/length);
00544
00545
00546 phioff[0] = phioffset + stphi[0];
00547 phioff[1] = phioffset + stphi[1];
00548
00549 for(int kk = 0; kk<nCell ; kk++){
00550
00551 inlow = phioff[0] + phi_bin*kk - phi_bin*0.5;
00552 inlow = fmod( inlow + 4.0 * pi , 2.0 * pi);
00553 inup = phioff[0] + phi_bin*kk + phi_bin*0.5;
00554 inup = fmod( inup + 4.0 * pi , 2.0 * pi);
00555 outlow = phioff[1] + phi_bin*kk - phi_bin*0.5;
00556 outlow = fmod( outlow + 4.0 * pi ,2.0 * pi);
00557 outup = phioff[1] + phi_bin*kk + phi_bin*0.5;
00558 outup = fmod( outup + 4.0 * pi , 2.0 * pi);
00559
00560 if((phi_in>=inlow && phi_in<=inup)) cellId_in = kk;
00561 if((phi_out>=outlow&&phi_out<outup)) cellId_out = kk;
00562 if(inlow > inup ){
00563 if((phi_in>=inlow&&phi_in<2.0*pi)||(phi_in>=0.0&&phi_in<inup)) cellId_in = kk;
00564 }
00565 if(outlow>outup){
00566 if((phi_out>=outlow&&phi_out<=2.0*pi)||(phi_out>=0.0&&phi_out<outup)) cellId_out = kk;
00567 }
00568 }
00569
00570 return (cellId_in==cellId_out);
00571 }
00572
00573 HepPoint3D MdcUtilitySvc::Hel(HepPoint3D piv, double dr,double phi0,double Alpha_L,double kappa,double dz,double dphi,double tanl) const{
00574 double x = piv.x() + dr*cos(phi0) + (Alpha_L/kappa) * (cos(phi0) - cos(phi0+dphi));
00575 double y = piv.y() + dr*sin(phi0) + (Alpha_L/kappa) * (sin(phi0) - sin(phi0+dphi));
00576 double z = piv.z() + dz - (Alpha_L/kappa) * dphi * tanl;
00577
00578 if((x>-1000. && x<1000.) || (y >-1000. && y <1000. ) ||(z>-1000. && z<1000.)){
00579 return HepPoint3D(x, y, z);
00580 }else{
00581 return HepPoint3D(0,0,0);
00582 }
00583 }
00584
00585
00586
00587 HepVector MdcUtilitySvc::besPar2PatPar(const HepVector& helixPar) const{
00588 HepVector helix(5,0);
00589 double d0 = -helixPar[0];
00590 double phi0 = helixPar[1]+ CLHEP::halfpi;
00591 double omega = Bz()*helixPar[2]/-333.567;
00592 double z0 = helixPar[3];
00593 double tanl = helixPar[4];
00594 helix[0] = d0;
00595 helix[1] = phi0;
00596 helix[2] = omega;
00597 helix[3] = z0;
00598 helix[4] = tanl;
00599 return helix;
00600 }
00601
00602
00603 HepSymMatrix MdcUtilitySvc::besErr2PatErr(const HepSymMatrix& err) const{
00604
00605
00606
00607
00608 HepSymMatrix mS(err.num_row(),0);
00609 mS[0][0]=-1.;
00610 mS[1][1]=1.;
00611 mS[2][2]=Bz()/-333.567;
00612 mS[3][3]=1.;
00613 mS[4][4]=1.;
00614 HepSymMatrix mVy= err.similarity(mS);
00615
00616 return mVy;
00617 }
00618
00619 double MdcUtilitySvc::p_cms(HepVector helix, int runNo, double mass) const{
00620 HepLorentzVector p4;
00621 p4.setPx(- sin(helix[1]) / fabs(helix[2]));
00622 p4.setPy(cos(helix[1]) / fabs(helix[2]));
00623 p4.setPz(helix[4] / fabs(helix[2]));
00624 double p3 = p4.mag();
00625 mass = 0.000511;
00626 p4.setE(sqrt(p3 * p3 + mass * mass));
00627
00628 double ecm;
00629 if (runNo > 9815) {
00630 ecm = 3.097;
00631 }else{
00632 ecm = 3.68632;
00633 }
00634 double zboost = 0.0075;
00635
00636 HepLorentzVector psip(0.011 * ecm, 0, zboost, ecm);
00637
00638 Hep3Vector boostv = psip.boostVector();
00639
00640
00641 p4.boost(- boostv);
00642
00643
00644 double p_cms = p4.rho();
00645
00646
00647
00648 return p_cms;
00649 }
00650
00651 Hep3Vector MdcUtilitySvc::momentum(const RecMdcTrack* trk) const{
00652
00653 double fi0 = trk->helix(1);
00654 double cpa = trk->helix(2);
00655 double tanl = trk->helix(4);
00656
00657 double pxy = 0.;
00658 if(cpa != 0) pxy = 1/fabs(cpa);
00659
00660 double px = pxy * (-sin(fi0));
00661 double py = pxy * cos(fi0);
00662 double pz = pxy * tanl;
00663
00664 Hep3Vector p;
00665 p.set(px, py, pz);
00666 return p;
00667 }
00668
00669 double MdcUtilitySvc::probab(const int& ndof, const double& chisq) const{
00670
00671
00672 double srtopi=2.0/sqrt(2.0*M_PI);
00673 double upl=100.0;
00674
00675 double prob=0.0;
00676 if(ndof<=0) {return prob;}
00677 if(chisq<0.0) {return prob;}
00678 if(ndof<=60) {
00679
00680 if(chisq>upl) {return prob;}
00681 double sum=exp(-0.5*chisq);
00682 double term=sum;
00683
00684 int m=ndof/2;
00685 if(2*m==ndof){
00686 if(m==1){return sum;}
00687 for(int i=2; i<=m;i++){
00688 term=0.5*term*chisq/(i-1);
00689 sum+=term;
00690 }
00691 return sum;
00692
00693
00694 }else{
00695
00696 double srty=sqrt(chisq);
00697 double temp=srty/M_SQRT2;
00698 prob=erfc(temp);
00699 if(ndof==1) {return prob;}
00700 if(ndof==3) {return (srtopi*srty*sum+prob);}
00701 m=m-1;
00702 for(int i=1; i<=m; i++){
00703 term=term*chisq/(2*i+1);
00704 sum+=term;
00705 }
00706 return (srtopi*srty*sum+prob);
00707
00708 }
00709
00710 }else{
00711
00712 double srty=sqrt(chisq)-sqrt(ndof-0.5);
00713 if(srty<12.0) {prob=0.5*erfc(srty);};
00714
00715 }
00716
00717 return prob;
00718 }
00719