MdcUtilitySvc Class Reference

#include <MdcUtilitySvc.h>

Inheritance diagram for MdcUtilitySvc:

IMdcUtilitySvc List of all members.

Public Member Functions

 MdcUtilitySvc (const std::string &name, ISvcLocator *svcloc)
 ~MdcUtilitySvc ()
virtual StatusCode initialize ()
virtual StatusCode finalize ()
StatusCode queryInterface (const InterfaceID &riid, void **ppvUnknown)
int nLayerTrackPassed (const HepVector helix) const
int nLayerTrackPassed (const double helix[5]) const
HepVector patPar2BesPar (const HepVector &helixPar) const
HepSymMatrix patErr2BesErr (const HepSymMatrix &err) const
HepVector besPar2PatPar (const HepVector &helixPar) const
HepSymMatrix besErr2PatErr (const HepSymMatrix &err) const
double doca (int layer, int cell, const HepVector helix, const HepSymMatrix errMat, bool passCellRequired=true, bool doSag=true) const
double doca (int layer, int cell, HepPoint3D eastP, HepPoint3D westP, const HepVector helixBes, const HepSymMatrix errMatBes, bool passCellRequired=true, bool doSag=true) const
double doca (int layer, int cell, const MdcSWire *sWire, const HepVector helixPat, const HepSymMatrix errMatPat, bool passCellRequired=true) const
double docaPatPar (int layer, int cell, const HepVector helixPat, const HepSymMatrix errMatPat, bool passCellRequired=true, bool doSag=true) const
double docaPatPar (int layer, int cell, HepPoint3D eastP, HepPoint3D westP, const HepVector helixBes, const HepSymMatrix errMatBes, bool passCellRequired=true, bool doSag=true) const
double docaPatPar (int layer, int cell, const MdcSWire *sWire, const HepVector helixPat, const HepSymMatrix errMatPat, bool passCellRequired=true) const
HepPoint3D pointOnHelix (const HepVector helixPar, int lay, int innerOrOuter) const
HepPoint3D pointOnHelixPatPar (const HepVector helixPat, int lay, int innerOrOuter) const
bool cellTrackPassedByPhi (const HepVector helix, int layer, int &cellId_in, int &cellId_out) const
bool cellTrackPassedByPhiPatPar (const HepVector helix, int layer, int &cellId_in, int &cellId_out) const
bool cellTrackPassed (const HepVector helix, int layer, int &cellId_in, int &cellId_out) const
bool cellTrackPassedPatPar (const HepVector helix, int layer, int &cellId_in, int &cellId_out) const
HepPoint3D Hel (HepPoint3D piv, double dr, double phi0, double Alpha_L, double kappa, double dz, double dphi, double tanl) const
double p_cms (HepVector helix, int runNo, double mass) const
Hep3Vector momentum (const RecMdcTrack *trk) const
double probab (const int &ndof, const double &chisq) const

Static Public Member Functions

static const InterfaceID & interfaceID ()

Private Member Functions

double Bz () const

Private Attributes

IMdcGeomSvc * m_mdcGeomSvc
IMagneticFieldSvc * m_pIMF
int m_debug
bool m_doSag

Detailed Description

Definition at line 19 of file MdcUtilitySvc.h.

Constructor & Destructor Documentation

MdcUtilitySvc::MdcUtilitySvc ( const std::string &  name,
ISvcLocator *  svcloc 
)

MdcUtilitySvc::~MdcUtilitySvc (  ) 

Definition at line 27 of file MdcUtilitySvc.cxx.

00027                              {
00028 }

Member Function Documentation

HepSymMatrix MdcUtilitySvc::besErr2PatErr ( const HepSymMatrix &  err  )  const [virtual]

Implements IMdcUtilitySvc.

Definition at line 603 of file MdcUtilitySvc.cxx.

References Bz().

Referenced by doca().

00603                                                                       {
00604   //error matrix
00605   //std::cout<<" err1  "<<err<<" "<<err.num_row()<<std::endl;
00606   //V(Y)=S * V(X) * ST , mS = S , mVy = V(Y) , err() = V(X)
00607   //int n = err.num_row();
00608   HepSymMatrix mS(err.num_row(),0);
00609   mS[0][0]=-1.;//dr0=-d0
00610   mS[1][1]=1.;
00611   mS[2][2]=Bz()/-333.567;//pxy -> cpa
00612   mS[3][3]=1.;
00613   mS[4][4]=1.;
00614   HepSymMatrix mVy= err.similarity(mS);
00615   //std::cout<<" err2  "<<n<<" "<<mVy<<std::endl;
00616   return mVy;
00617 }

HepVector MdcUtilitySvc::besPar2PatPar ( const HepVector &  helixPar  )  const [virtual]

Implements IMdcUtilitySvc.

Definition at line 587 of file MdcUtilitySvc.cxx.

References Bz(), and phi0.

Referenced by cellTrackPassedByPhi(), doca(), MdcxCosmicSewer::execute(), nLayerTrackPassed(), and pointOnHelix().

00587                                                                      {
00588   HepVector helix(5,0);
00589   double d0 = -helixPar[0];    //cm
00590   double phi0 = helixPar[1]+ CLHEP::halfpi;
00591   double omega = Bz()*helixPar[2]/-333.567;
00592   double z0 = helixPar[3];    //cm
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 }

double MdcUtilitySvc::Bz (  )  const [inline, private]

Definition at line 57 of file MdcUtilitySvc.h.

Referenced by besErr2PatErr(), besPar2PatPar(), cellTrackPassed(), nLayerTrackPassed(), patErr2BesErr(), and patPar2BesPar().

00057 { return m_pIMF->getReferField()*1000.; };

bool MdcUtilitySvc::cellTrackPassed ( const HepVector  helix,
int  layer,
int &  cellId_in,
int &  cellId_out 
) const [virtual]

Implements IMdcUtilitySvc.

Definition at line 425 of file MdcUtilitySvc.cxx.

References Bz(), cos(), Hel(), IMdcGeomSvc::Layer(), MdcGeoLayer::Length(), m_mdcGeomSvc, MdcGeoLayer::NCell(), MdcGeoLayer::Offset(), phi0, Constants::pi, pi, MdcGeoLayer::Radius(), MdcGeoLayer::RCSiz1(), MdcGeoLayer::RCSiz2(), MdcGeoLayer::Shift(), sign(), MdcGeoLayer::Slant(), MdcGeoLayer::Sup(), MdcGeoSuper::Type(), type, and x.

Referenced by cellTrackPassedPatPar().

00425                                                                                                    {
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()); //magnetic field const
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   //retrieve Mdc geometry information
00474   rCize1 = 0.1 * m_mdcGeomSvc->Layer(lay)->RCSiz1();  //mm -> cm 
00475   rCize2 = 0.1 * m_mdcGeomSvc->Layer(lay)->RCSiz2();  //mm -> cm
00476   rLay   = 0.1 * m_mdcGeomSvc->Layer(lay)->Radius();  //mm -> cm
00477   length = 0.1 * m_mdcGeomSvc->Layer(lay)->Length();  //mm -> cm
00478   //double halfLength=length/2.;
00479   //(conversion of the units mm(mc)->cm(rec))
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 ; //radius of inner field wir ( beam pipe)
00487   m_crio[1] = rLay + rCize2 ; //radius of outer field wir ( beam pipe)
00488 
00489   int sign = -1; //assumed the first half circle
00490 
00491   Hep3Vector iocand[2];
00492   Hep3Vector cell_IO[2];
00493 
00494   for(int ii =0; ii<2; ii++){
00495     // By law of cosines to calculate the alpha(up and down field wir_Ge)
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] ){  //very very nausea
00524       dphi[ii] -= 2*pi;
00525     }
00526 
00527     cell_IO[ii] = Hel(piv, dr0, phi, ALPHA_loc, kappa,dz0, dphi[ii], tanl); 
00528 
00529     //cout<<" cell_IO["<<ii<<"] : "<<cell_IO[ii]<<endl;
00530     if( (cell_IO[ii].x()==0 ) && (cell_IO[ii].y()==0) && (cell_IO[ii].z()==0)) continue;
00531   }
00532   //if(((fabs(cell_IO[0].z())*10-halfLength)>-7.) && ((fabs(cell_IO[1].z())*10-halfLength)>-7.))return true; //Out sensitive area
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   //determine the in/out cell id
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   //stphi[0],stphi[1] to position fo z axis ,the angle of deflxsion in x_Y
00545 
00546   phioff[0] = phioffset + stphi[0];
00547   phioff[1] = phioffset + stphi[1];
00548 
00549   for(int kk = 0; kk<nCell ; kk++){
00550     //for stereo layer
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   }//end of nCell loop  
00569 
00570   return (cellId_in==cellId_out);
00571 }

bool MdcUtilitySvc::cellTrackPassedByPhi ( const HepVector  helix,
int  layer,
int &  cellId_in,
int &  cellId_out 
) const [virtual]

Implements IMdcUtilitySvc.

Definition at line 378 of file MdcUtilitySvc.cxx.

References besPar2PatPar(), and cellTrackPassedByPhiPatPar().

00378                                                                                                                {
00379   HepVector helixPat = besPar2PatPar(helixBes);
00380   return cellTrackPassedByPhiPatPar(helixPat, layer, cellId_in, cellId_out);
00381 }

bool MdcUtilitySvc::cellTrackPassedByPhiPatPar ( const HepVector  helix,
int  layer,
int &  cellId_in,
int &  cellId_out 
) const [virtual]

Implements IMdcUtilitySvc.

Definition at line 387 of file MdcUtilitySvc.cxx.

References MdcGeoWire::Backward(), cos(), MdcGeoWire::Forward(), IMdcGeomSvc::Layer(), m_debug, m_mdcGeomSvc, MdcGeoLayer::NCell(), BesAngle::rad(), and IMdcGeomSvc::Wire().

Referenced by cellTrackPassedByPhi(), and docaPatPar().

00387                                                                                                                      {
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.;//mm2cm 
00392   double rMid = rEnd * cos(dPhiz);
00393   //std::cout<<"( cell "<<0<<" westPphi "<<m_mdcGeomSvc->Wire(layer,0)->Forward().phi() <<" eastPphi "
00394   //<<m_mdcGeomSvc->Wire(layer,0)->Backward().phi()<<" twist  "<<dPhiz<<" rend "<<rEnd<<std::endl;
00395   double fltLenIn = rMid;
00396   double phiIn = helixPat[1] + helixPat[2]*fltLenIn;//phi0 + omega * L
00397 
00398   double phiEPOffset = m_mdcGeomSvc->Wire(layer,0)->Backward().phi();//east.phi()= 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   //BesAngle tmp(phiIn - layPtr->phiEPOffset());
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;//pass more than one cell
00415 }

bool MdcUtilitySvc::cellTrackPassedPatPar ( const HepVector  helix,
int  layer,
int &  cellId_in,
int &  cellId_out 
) const [virtual]

Implements IMdcUtilitySvc.

Definition at line 418 of file MdcUtilitySvc.cxx.

References cellTrackPassed(), and patPar2BesPar().

00418                                                                                                            {
00419   HepVector helixBes = patPar2BesPar(helixPat);
00420   return cellTrackPassed(helixBes, layer, cellId_in, cellId_out);
00421 }

double MdcUtilitySvc::doca ( int  layer,
int  cell,
const MdcSWire *  sWire,
const HepVector  helixPat,
const HepSymMatrix  errMatPat,
bool  passCellRequired = true 
) const [virtual]

Implements IMdcUtilitySvc.

Definition at line 234 of file MdcUtilitySvc.cxx.

References besErr2PatErr(), besPar2PatPar(), and docaPatPar().

00234                                                                                                                                                         {
00235   HepVector helixPat = besPar2PatPar(helixBes);
00236   HepSymMatrix errMatPat= besErr2PatErr(errMatBes);
00237 
00238   return docaPatPar(layer, cell, sWire, helixPat, errMatPat, passCellRequired);//call 6.
00239 }

double MdcUtilitySvc::doca ( int  layer,
int  cell,
HepPoint3D  eastP,
HepPoint3D  westP,
const HepVector  helixBes,
const HepSymMatrix  errMatBes,
bool  passCellRequired = true,
bool  doSag = true 
) const [virtual]

Implements IMdcUtilitySvc.

Definition at line 225 of file MdcUtilitySvc.cxx.

References besErr2PatErr(), besPar2PatPar(), and docaPatPar().

00225                                                                                                                                                                                 {
00226   HepVector helixPat = besPar2PatPar(helixBes);
00227   HepSymMatrix errMatPat= besErr2PatErr(errMatBes);
00228 
00229   return docaPatPar(layer, cell, eastP, westP, helixPat, errMatPat, passCellRequired, doSag);//call 5.
00230 }

double MdcUtilitySvc::doca ( int  layer,
int  cell,
const HepVector  helix,
const HepSymMatrix  errMat,
bool  passCellRequired = true,
bool  doSag = true 
) const [virtual]

Implements IMdcUtilitySvc.

Definition at line 216 of file MdcUtilitySvc.cxx.

References besErr2PatErr(), besPar2PatPar(), and docaPatPar().

Referenced by docaPatPar(), EsTimeAlg::execute(), MilleAlign::fillHist(), MilleAlign::getDeriGlo(), and MilleAlign::getDeriLoc().

00216                                                                                                                                               {
00217   HepVector helixPat = besPar2PatPar(helixBes);
00218   HepSymMatrix errMatPat= besErr2PatErr(errMatBes);
00219 
00220   return docaPatPar(layer, cell, helixPat, errMatPat, passCellRequired, doSag);//call 4.
00221 }

double MdcUtilitySvc::docaPatPar ( int  layer,
int  cell,
const MdcSWire *  sWire,
const HepVector  helixPat,
const HepSymMatrix  errMatPat,
bool  passCellRequired = true 
) const [virtual]

Implements IMdcUtilitySvc.

Definition at line 277 of file MdcUtilitySvc.cxx.

References cellTrackPassedByPhiPatPar(), TrkPoca::doca(), doca(), TrkPocaBase::flt2(), MdcSWire::getEastPoint(), MdcSWire::getTraj(), Trajectory::hiRange(), Trajectory::lowRange(), m_debug, m_mdcGeomSvc, pointOnHelixPatPar(), MdcSagTraj::sag(), MdcSWire::xWireDC(), MdcSWire::yWireDC(), and MdcSWire::zEndDC().

00277                                                                                                                                                            {
00278 
00279   if(m_debug) std::cout<<" helixPat  "<<helixPat<<std::endl;
00280   if(m_debug) std::cout<<" layer "<<layer<<" cell "<<cell<<std::endl;
00281   //-----cell track passed
00282   int cellId_in  = -1;
00283   int cellId_out = -1;
00284 
00285   //cellTrackPassedPatPar(helixPat,layer,cellId_in,cellId_out);//yzhang FIXME 2012-07-11 
00286   cellTrackPassedByPhiPatPar(helixPat,layer,cellId_in,cellId_out);//yzhang FIXME 2012-07-11 
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   //-----helix trajectory
00295   HelixTraj* helixTraj = new HelixTraj(helixPat,errMat);
00296   const Trajectory* trajHelix    = dynamic_cast<const Trajectory*> (helixTraj);
00297 
00298   //-----pointOnHelix
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   //------wire trajectory
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   //const HepPoint3D* stop_In   = sWire->getWestPoint();
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   //std::cout<< " west -------- "<< stop_In->x()<<","<<stop_In->y()<<","<<stop_In->z() <<std::endl;
00317 
00318   //------calc poca
00319   double doca  = -999.;
00320   TrkPoca* trkPoca;
00321   double zWire = cell_IO.z();
00322   //HepPoint3D pos_in(zWire,sWire->xWireDC(zWire),sWire->yWireDC(zWire)) ;
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   //if(m_debug) std::cout<<"pos_in "<<pos_in<<" fltWire  "<<fltWire<<std::endl;
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   //std::cout<<" hitLen  "<<hitLen<<" startLen  "<<startLen<<" endLen  "<<endLen <<" doca "<<doca<<std::endl;
00341 
00342   if(m_debug) std::cout<<" doca  "<<doca<<std::endl;
00343   delete helixTraj;
00344   delete trkPoca;
00345 
00346   return doca;
00347 
00348 } //----------end of calculatng the doca ---------------------------------//

double MdcUtilitySvc::docaPatPar ( int  layer,
int  cell,
HepPoint3D  eastP,
HepPoint3D  westP,
const HepVector  helixBes,
const HepSymMatrix  errMatBes,
bool  passCellRequired = true,
bool  doSag = true 
) const [virtual]

Implements IMdcUtilitySvc.

Definition at line 261 of file MdcUtilitySvc.cxx.

References doca(), docaPatPar(), MdcGeoWire::Id(), m_mdcGeomSvc, MdcGeoWire::Sag(), and IMdcGeomSvc::Wire().

00261                                                                                                                                                                                   {
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.;//mm2cm
00267   const MdcSWire* sWire = new MdcSWire(eastP, westP, sag, id , cell);//cm
00268 
00269   double doca = docaPatPar(layer, cell, sWire, helixPat, errMat, passCellRequired);//call 6.
00270 
00271   delete sWire;
00272   return doca;
00273 }

double MdcUtilitySvc::docaPatPar ( int  layer,
int  cell,
const HepVector  helixPat,
const HepSymMatrix  errMatPat,
bool  passCellRequired = true,
bool  doSag = true 
) const [virtual]

Implements IMdcUtilitySvc.

Definition at line 243 of file MdcUtilitySvc.cxx.

References MdcGeoWire::Backward(), doca(), MdcGeoWire::Forward(), MdcGeoWire::Id(), m_mdcGeomSvc, MdcGeoWire::Sag(), and IMdcGeomSvc::Wire().

Referenced by doca(), docaPatPar(), and MdcxCosmicSewer::execute().

00243                                                                                                                                                  {
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.;//mm2cm
00249   HepPoint3D eastP = geowir->Backward()/10.0;//mm2cm
00250   HepPoint3D westP = geowir->Forward() /10.0;//mm2cm
00251   const MdcSWire* sWire = new MdcSWire(eastP, westP, sag, id , cell);
00252 
00253   double doca = docaPatPar(layer, cell, sWire, helixPat, errMat, passCellRequired);//call 6.
00254 
00255   delete sWire;
00256   return doca;
00257 }

StatusCode MdcUtilitySvc::finalize (  )  [virtual]

Definition at line 57 of file MdcUtilitySvc.cxx.

References Bes_Common::INFO.

00057                                   {
00058   MsgStream log(messageService(), name());
00059   log << MSG::INFO << "MdcUtilitySvc::finalize()" << endreq;
00060 
00061   return StatusCode::SUCCESS;
00062 }

HepPoint3D MdcUtilitySvc::Hel ( HepPoint3D  piv,
double  dr,
double  phi0,
double  Alpha_L,
double  kappa,
double  dz,
double  dphi,
double  tanl 
) const [virtual]

Implements IMdcUtilitySvc.

Definition at line 573 of file MdcUtilitySvc.cxx.

References cos(), sin(), and x.

Referenced by cellTrackPassed().

00573                                                                                                                                       {
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   //cout<<"HepPoint3D(x, y, z) = "<<HepPoint3D(x, y, z)<<endl;
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 }

StatusCode MdcUtilitySvc::initialize (  )  [virtual]

Definition at line 30 of file MdcUtilitySvc.cxx.

References calibUtil::ERROR, Bes_Common::FATAL, Bes_Common::INFO, m_mdcGeomSvc, and m_pIMF.

00030                                     {
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   //Initalze magnetic flied 
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   // Initialize MdcGeomSvc
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 }

static const InterfaceID& IMdcUtilitySvc::interfaceID (  )  [inline, static, inherited]

Definition at line 18 of file IMdcUtilitySvc.h.

References IID_IMdcUtilitySvc().

00018 { return IID_IMdcUtilitySvc; }

Hep3Vector MdcUtilitySvc::momentum ( const RecMdcTrack *  trk  )  const [virtual]

Implements IMdcUtilitySvc.

Definition at line 651 of file MdcUtilitySvc.cxx.

References cos(), DstMdcTrack::helix(), and sin().

00651                                                               {
00652   // double dr   = trk->helix(0);
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); // MeV/c
00666   return p;
00667 }

int MdcUtilitySvc::nLayerTrackPassed ( const double  helix[5]  )  const [virtual]

Implements IMdcUtilitySvc.

Definition at line 77 of file MdcUtilitySvc.cxx.

References genRecEmupikp::i, and nLayerTrackPassed().

00077                                                                   {
00078 
00079   HepVector helixBesParam(5,0);
00080   for(int i=0; i<5; ++i) helixBesParam[i] = helixBes[i];
00081 
00082   return nLayerTrackPassed(helixBesParam);
00083 }

int MdcUtilitySvc::nLayerTrackPassed ( const HepVector  helix  )  const [virtual]

Implements IMdcUtilitySvc.

Definition at line 86 of file MdcUtilitySvc.cxx.

References alpha, besPar2PatPar(), Bz(), IMdcGeomSvc::Layer(), MdcGeoLayer::Length(), m_mdcGeomSvc, phi0, and MdcGeoLayer::Radius().

Referenced by nLayerTrackPassed().

00086                                                                   {
00087   HepVector helix= besPar2PatPar(helixBes);
00088 
00089   int nLayer = 0;
00090 
00091   for(unsigned iLayer=0; iLayer<43; iLayer++){
00092     //flightLength is the path length of track in the x-y plane
00093     //guess flightLength by the radius in middle of the wire.
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.; //unit from m/s to cm/s
00100     double alpha = 1/(c * Bz());//~333.567
00101     double kappa = helix[2];
00102     double rc = (-1.)*alpha/kappa;
00103     //std::cout<<"MdcUtilitySvc alpha   "<<alpha<<std::endl;
00104     double tanl = helix[4];
00105     double phi0 = helix[1];
00106     double phi = flightLength/rc + phi0;//turning angle
00107     double z = pivot.z() + dz - (alpha/kappa) * tanl * phi;
00108 
00109     double layerHalfLength = m_mdcGeomSvc->Layer(iLayer)->Length()/2.;
00110 
00111     //std::cout<<"MdcUtilitySvc length  "<<layerHalfLength<<std::endl;
00112 
00113     if (fabs(z) < fabs(layerHalfLength)) ++nLayer;
00114   }
00115 
00116   return nLayer;
00117 }

double MdcUtilitySvc::p_cms ( HepVector  helix,
int  runNo,
double  mass 
) const [virtual]

Implements IMdcUtilitySvc.

Definition at line 619 of file MdcUtilitySvc.cxx.

References cos(), and sin().

00619                                                                         {
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   //HepLorentzVector psip(0.011 * 3.68632, 0, 0.0075, 3.68632);
00636   HepLorentzVector psip(0.011 * ecm, 0, zboost, ecm);
00637   //cout << setw(15) << ecm << setw(15) << zboost << endl;
00638   Hep3Vector boostv = psip.boostVector();
00639 
00640   //std::cout<<__FILE__<<" boostv "<<boostv<<  std::endl;
00641   p4.boost(- boostv);
00642 
00643   //std::cout<<__FILE__<<" p4 "<<p4<<  std::endl;
00644   double p_cms = p4.rho();
00645   //phicms = p4.phi();
00646   //p4.boost(boostv);
00647 
00648   return p_cms;
00649 }

HepSymMatrix MdcUtilitySvc::patErr2BesErr ( const HepSymMatrix &  err  )  const [virtual]

Implements IMdcUtilitySvc.

Definition at line 137 of file MdcUtilitySvc.cxx.

References Bz().

00137                                                                       {
00138   //error matrix
00139   //std::cout<<" err1  "<<err<<" "<<err.num_row()<<std::endl;
00140   //V(Y)=S * V(X) * ST , mS = S , mVy = V(Y) , err() = V(X)
00141   //int n = err.num_row();
00142   HepSymMatrix mS(err.num_row(),0);
00143   mS[0][0]=-1.;//dr0=-d0
00144   mS[1][1]=1.;
00145 
00146   mS[2][2]=Bz()/-333.567;//pxy -> cpa
00147   mS[3][3]=1.;
00148   mS[4][4]=1.;
00149   HepSymMatrix mVy= err.similarity(mS);
00150   //std::cout<<" err2  "<<n<<" "<<mVy<<std::endl;
00151   return mVy;
00152 }

HepVector MdcUtilitySvc::patPar2BesPar ( const HepVector &  helixPar  )  const [virtual]

Implements IMdcUtilitySvc.

Definition at line 120 of file MdcUtilitySvc.cxx.

References Bz(), and phi0.

Referenced by cellTrackPassedPatPar().

00120                                                                      {
00121   HepVector helix(5,0);
00122   double d0 = -helixPar[0];    //cm
00123   double phi0 = helixPar[1]+ CLHEP::halfpi;
00124   double omega = Bz()*helixPar[2]/-333.567;
00125   double z0 = helixPar[3];    //cm
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   //std::cout<<"helix   "<<helix<<std::endl;
00133   return helix;
00134 }

HepPoint3D MdcUtilitySvc::pointOnHelix ( const HepVector  helixPar,
int  lay,
int  innerOrOuter 
) const [virtual]

Implements IMdcUtilitySvc.

Definition at line 156 of file MdcUtilitySvc.cxx.

References besPar2PatPar(), and pointOnHelixPatPar().

00156                                                                                                  {
00157   HepVector helixPat= besPar2PatPar(helixBes);
00158   return pointOnHelixPatPar(helixPat, layer, innerOrOuter);
00159 }

HepPoint3D MdcUtilitySvc::pointOnHelixPatPar ( const HepVector  helixPat,
int  lay,
int  innerOrOuter 
) const [virtual]

Implements IMdcUtilitySvc.

Definition at line 163 of file MdcUtilitySvc.cxx.

References abs, cos(), Constants::epsilon, IMdcGeomSvc::Layer(), m_debug, m_mdcGeomSvc, phi0, MdcGeoLayer::Radius(), MdcGeoLayer::RCSiz1(), MdcGeoLayer::RCSiz2(), sin(), and x.

Referenced by docaPatPar(), and pointOnHelix().

00163                                                                                                        {
00164 
00165   double rInner,rOuter;
00166   //retrieve Mdc geometry information
00167   double rCize1 =0.1* m_mdcGeomSvc->Layer(layer)->RCSiz1();  //mm -> cm 
00168   double rCize2 =0.1* m_mdcGeomSvc->Layer(layer)->RCSiz2();  //mm -> cm
00169   double rLay   =0.1* m_mdcGeomSvc->Layer(layer)->Radius();  //mm -> cm
00170 
00171   //(conversion of the units mm(mc)->cm(rec))
00172   rInner = rLay - rCize1 ; //radius of inner field wire
00173   rOuter = rLay + rCize2 ; //radius of outer field wire
00174 
00175   //int sign = -1; //assumed the first half circle
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   //std::cout<<" xc yc "<<xc<<" "<<yc
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 }

double MdcUtilitySvc::probab ( const int &  ndof,
const double &  chisq 
) const [virtual]

Implements IMdcUtilitySvc.

Definition at line 669 of file MdcUtilitySvc.cxx.

References exp(), genRecEmupikp::i, M_PI, and subSeperate::temp.

00669                                                                       {
00670 
00671   //constants
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     //full treatment
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       }//(int i=2; i<=m)
00691       return sum;
00692       //even
00693 
00694     }else{
00695       //odd
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       }//(int i=1; i<=m; i++)
00706       return (srtopi*srty*sum+prob);
00707 
00708     }//(2*m==ndof)
00709 
00710   }else{
00711     //asymtotic Gaussian approx
00712     double srty=sqrt(chisq)-sqrt(ndof-0.5);
00713     if(srty<12.0) {prob=0.5*erfc(srty);};
00714 
00715   }//ndof<30
00716 
00717   return prob;
00718 }//endof probab

StatusCode MdcUtilitySvc::queryInterface ( const InterfaceID &  riid,
void **  ppvUnknown 
)

Definition at line 64 of file MdcUtilitySvc.cxx.

References IID_IMdcUtilitySvc().

00064                                                                                     {
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 }

Member Data Documentation

int MdcUtilitySvc::m_debug [private]

Definition at line 61 of file MdcUtilitySvc.h.

Referenced by cellTrackPassedByPhiPatPar(), docaPatPar(), and pointOnHelixPatPar().

bool MdcUtilitySvc::m_doSag [private]

Definition at line 62 of file MdcUtilitySvc.h.

IMdcGeomSvc* MdcUtilitySvc::m_mdcGeomSvc [private]

Definition at line 57 of file MdcUtilitySvc.h.

Referenced by cellTrackPassed(), cellTrackPassedByPhiPatPar(), docaPatPar(), initialize(), nLayerTrackPassed(), and pointOnHelixPatPar().

IMagneticFieldSvc* MdcUtilitySvc::m_pIMF [private]

Definition at line 60 of file MdcUtilitySvc.h.

Referenced by initialize().

Generated on Tue Nov 29 23:20:20 2016 for BOSS_7.0.2 by  doxygen 1.4.7