Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

MdcCheckUtil Class Reference

#include <MdcCheckUtil.h>

List of all members.

Public Member Functions

HepSymMatrix besErr2PatErr (const HepSymMatrix &err)
HepSymMatrix besErr2PatErr (const HepSymMatrix &err)
HepVector besPar2PatPar (const double helix[5])
HepVector besPar2PatPar (const HepVector &helixPar)
HepVector besPar2PatPar (const double helix[5])
HepVector besPar2PatPar (const HepVector &helixPar)
int cellOfR (HepVector helixBabar, int layer, double r)
int cellOfR (HepVector helixBabar, int layer, double r)
void cellPassed (const RecMdcTrack *tk, bool tkPass[43][288])
void cellPassed (const RecMdcTrack *tk, bool tkPass[43][288])
int cellTrackPassed (const HepVector helix, int layer, int &cellId_in, int &cellId_out)
int cellTrackPassed (const HepVector helix, int layer, int &cellId_in, int &cellId_out)
bool cellTrackPassedBelle (HepVector helix, int lay, int &cellId_in, int &cellId_cout)
bool cellTrackPassedBelle (HepVector helix, int lay, int &cellId_in, int &cellId_cout)
double doca (int layer, int cell, const MdcSWire *sWire, const HepVector helix, const HepSymMatrix errMat)
double doca (int layer, int cell, HepPoint3D eastP, HepPoint3D westP, const HepVector helix, const HepSymMatrix errMat)
double doca (int layer, int cell, const HepVector helix, const HepSymMatrix errMat)
double doca (int layer, int cell, const MdcSWire *sWire, const HepVector helix, const HepSymMatrix errMat)
double doca (int layer, int cell, HepPoint3D eastP, HepPoint3D westP, const HepVector helix, const HepSymMatrix errMat)
double doca (int layer, int cell, const HepVector helix, const HepSymMatrix errMat)
double docaPatPar (int layer, int cell, const MdcSWire *sWire, const HepVector helixPat, const HepSymMatrix errMatPat)
double docaPatPar (int layer, int cell, const HepVector helixPat, const HepSymMatrix errMatPat)
double docaPatPar (int layer, int cell, const MdcSWire *sWire, const HepVector helixPat, const HepSymMatrix errMatPat)
double docaPatPar (int layer, int cell, const HepVector helixPat, const HepSymMatrix errMatPat)
void dumpRecMdcTrack ()
void dumpRecMdcTrack ()
HepPoint3D Hel (HepPoint3D piv, double dr, double phi0, double Alpha_L, double kappa, double dz, double dphi, double tanl)
HepPoint3D Hel (HepPoint3D piv, double dr, double phi0, double Alpha_L, double kappa, double dz, double dphi, double tanl)
 MdcCheckUtil (bool doSag=false)
 MdcCheckUtil (bool doSag=false)
Hep3Vector momentum (const RecMdcTrack *trk)
Hep3Vector momentum (const RecMdcTrack *trk)
int nLayerTrackPassed (const double helix[5])
int nLayerTrackPassed (const HepVector helix)
int nLayerTrackPassed (const double helix[5])
int nLayerTrackPassed (const HepVector helix)
double p_cms (HepVector helix, int runNo, double mass)
double p_cms (HepVector helix, int runNo, double mass)
HepSymMatrix patRecErr2BesErr (const HepSymMatrix &err)
HepSymMatrix patRecErr2BesErr (const HepSymMatrix &err)
HepVector patRecPar2BesPar (const HepVector &helixPar)
HepVector patRecPar2BesPar (const HepVector &helixPar)
HepPoint3D pointOnHelix (const HepVector helixPar, int lay, int innerOrOuter)
HepPoint3D pointOnHelix (const HepVector helixPar, int lay, int innerOrOuter)
HepPoint3D pointOnHelixPatPar (const HepVector helixPat, int lay, int innerOrOuter)
HepPoint3D pointOnHelixPatPar (const HepVector helixPat, int lay, int innerOrOuter)
double probab (const int &ndof, const double &chisq)
double probab (const int &ndof, const double &chisq)
void setPassCellRequired (bool pass)
void setPassCellRequired (bool pass)
 ~MdcCheckUtil ()
 ~MdcCheckUtil ()

Private Attributes

double Bz
int m_debug
const MdcDetector * m_gm
const MdcDetector * m_gm
MdcGeomSvc * m_mdcGeomSvc
MdcGeomSvc * m_mdcGeomSvc
bool m_passCellRequired
IMagneticFieldSvc * m_pIMF
IMagneticFieldSvc * m_pIMF

Constructor & Destructor Documentation

MdcCheckUtil::MdcCheckUtil (  bool  doSag = false  ) 
  

00040                                     {
00041 
00042   //Initalze magnetic flied 
00043   IService* svc;
00044   Gaudi::svcLocator()->getService("MagneticFieldSvc",svc);
00045   m_pIMF = dynamic_cast<IMagneticFieldSvc*> (svc);
00046   if(! m_pIMF){
00047     std::cout<<" ERROR Unable to open Magnetc feld service "<<std::endl;
00048   }
00049   //get Bz for Check TEMP, Bz may be changed by run
00050   Bz = m_pIMF->getReferField()*1000.;
00051 
00052   // Initialize MdcGeomSvc
00053   Gaudi::svcLocator()->getService("MdcGeomSvc",svc);
00054   m_mdcGeomSvc= dynamic_cast<MdcGeomSvc*> (svc);
00055   if(! m_mdcGeomSvc){
00056     std::cout<<" FATAL Could not load MdcGeomSvc!  "<<std::endl;
00057   }
00058 
00059 
00060   //Initialze MdcDetector
00061   m_gm = MdcDetector::instance(doSag);
00062 
00063   m_passCellRequired = true;
00064 }

MdcCheckUtil::~MdcCheckUtil (   )  [inline]
  

00018 {};

MdcCheckUtil::MdcCheckUtil (  bool  doSag = false  ) 
 

MdcCheckUtil::~MdcCheckUtil (   )  [inline]
  

00018 {};

Member Function Documentation

HepSymMatrix MdcCheckUtil::besErr2PatErr (  const HepSymMatrix &  err  ) 
 

HepSymMatrix MdcCheckUtil::besErr2PatErr (  const HepSymMatrix &  err  ) 
  

00531                                                                {
00532   //error matrix
00533   //std::cout<<" err1  "<<err<<" "<<err.num_row()<<std::endl;
00534   //V(Y)=S * V(X) * ST , mS = S , mVy = V(Y) , err() = V(X)
00535   //int n = err.num_row();
00536   HepSymMatrix mS(err.num_row(),0);
00537   mS[0][0]=-1.;//dr0=-d0
00538   mS[1][1]=1.;
00539   mS[2][2]=Bz/-333.567;//pxy -> cpa
00540   mS[3][3]=1.;
00541   mS[4][4]=1.;
00542   HepSymMatrix mVy= err.similarity(mS);
00543   //std::cout<<" err2  "<<n<<" "<<mVy<<std::endl;
00544   return mVy;
00545 }

HepVector MdcCheckUtil::besPar2PatPar (  const double  helix[5]  ) 
 

HepVector MdcCheckUtil::besPar2PatPar (  const HepVector &  helixPar  ) 
 

HepVector MdcCheckUtil::besPar2PatPar (  const double  helix[5]  ) 
  

00503                                                              {
00504   HepVector helix(5,0);
00505   helix[0]=helixPar[0];
00506   helix[1]=helixPar[1];
00507   helix[2]=helixPar[2];
00508   helix[3]=helixPar[3];
00509   helix[4]=helixPar[4];
00510   return besPar2PatPar(helix);
00511 }

HepVector MdcCheckUtil::besPar2PatPar (  const HepVector &  helixPar  ) 
  

00515                                                               {
00516   HepVector helix(5,0);
00517   double d0 = -helixPar[0];    //cm
00518   double phi0 = helixPar[1]+ CLHEP::halfpi;
00519   double omega = Bz*helixPar[2]/-333.567;
00520   double z0 = helixPar[3];    //cm
00521   double tanl = helixPar[4];
00522   helix[0] = d0;
00523   helix[1] = phi0;
00524   helix[2] = omega;
00525   helix[3] = z0;
00526   helix[4] = tanl;
00527   return helix;
00528 }

int MdcCheckUtil::cellOfR (  HepVector  helixBabar,
int  layer,
double  r
) 
 

int MdcCheckUtil::cellOfR (  HepVector  helixBabar,
int  layer,
double  r
) 
 

void MdcCheckUtil::cellPassed (  const RecMdcTrack *  tk,
bool  tkPass[43][288]
) 
 

void MdcCheckUtil::cellPassed (  const RecMdcTrack *  tk,
bool  tkPass[43][288]
) 
  

00314                                                                               {
00315   HepVector helix(5);
00316   helix[0]=tk->helix(0);
00317   helix[1]=tk->helix(1);
00318   helix[2]=tk->helix(2);
00319   helix[3]=tk->helix(3);
00320   helix[4]=tk->helix(4);
00321 
00322   for(int i=0;i<43;i++){
00323     for(int j=0; j<288; j++){
00324       cellTkPassed[i][j]=false;
00325     }
00326   }
00327 
00328   for (int layer=0; layer<43; layer++){
00329     //-----cell track passed
00330 
00331     int cellId_in  = -1;
00332     int cellId_out = -1;
00333     cellTrackPassedBelle(helix,layer,cellId_in,cellId_out);
00334     cellTkPassed[layer][cellId_in] = true;
00335     cellTkPassed[layer][cellId_out] = true;
00336   }
00337 
00338 }

int MdcCheckUtil::cellTrackPassed (  const HepVector  helix,
int  layer,
int &  cellId_in,
int &  cellId_out
) 
 

int MdcCheckUtil::cellTrackPassed (  const HepVector  helix,
int  layer,
int &  cellId_in,
int &  cellId_out
) 
  

00294                                                                                                {
00295   const MdcLayer* layPtr = m_gm->Layer(layer);
00296   int nCell = layPtr->nWires();
00297 
00298   double fltLenIn = layPtr->rMid(); 
00299   double phiIn = helix[1] + helix[2]*fltLenIn;//phi0 + omega * L
00300 
00301   BesAngle tmp(phiIn - layPtr->phiEPOffset());
00302   int wlow = (int)floor(nCell * tmp.rad() / twopi );
00303   int wbig = (int)ceil(nCell * tmp.rad() / twopi );
00304   if (tmp == 0 ){ wlow = -1; wbig = 1; }
00305   cellId_in  = mdcWrapWire(wlow, nCell);
00306   cellId_out = mdcWrapWire(wbig, nCell);
00307 
00308   bool passedOneCell = (cellId_in == cellId_out);
00309 
00310   return passedOneCell;//pass more than one cell
00311 }

bool MdcCheckUtil::cellTrackPassedBelle (  HepVector  helix,
int  lay,
int &  cellId_in,
int &  cellId_cout
) 
 

bool MdcCheckUtil::cellTrackPassedBelle (  HepVector  helix,
int  lay,
int &  cellId_in,
int &  cellId_cout
) 
  

00341                                                                                               {
00342   int charge,type,nCell;
00343   double dr0,phi0,kappa,dz0,tanl;
00344   double ALPHA_loc,rho,phi,cosphi0,sinphi0,x_hc,y_hc,z_hc;
00345   double dphi0,IO_phi,C_alpha,xx,yy;
00346   double inlow,inup,outlow,outup,phi_in,phi_out,phi_bin;
00347   double rCize1,rCize2,rLay,length,phioffset,slant,shift; 
00348   double m_crio[2], phi_io[2], stphi[2],phioff[2],dphi[2];
00349 
00350   dr0   = helix[0];  
00351   phi0  = helix[1];
00352   kappa = helix[2];
00353   dz0   = helix[3];
00354   tanl  = helix[4];
00355 
00356   ALPHA_loc = 1000/(2.99792458*Bz); //magnetic field const
00357   charge    = ( kappa >=0 ) ? 1 : -1 ;
00358   rho       = ALPHA_loc/kappa ;
00359   double pi = Constants::pi;
00360   phi      = fmod(phi0 + 4*pi , 2*pi);
00361   cosphi0  = cos(phi);
00362   sinphi0  = (1.0 - cosphi0 ) * (1.0 + cosphi0 );
00363   sinphi0  = sqrt(( sinphi0 > 0.) ? sinphi0 : 0.);
00364   if( phi > pi ) sinphi0 = -sinphi0 ;
00365 
00366   x_hc     = 0. + ( dr0 + rho ) * cosphi0;
00367   y_hc     = 0. + ( dr0 + rho ) * sinphi0;
00368   z_hc     = 0. + dz0;
00369 
00370 
00371   HepPoint3D piv(0.,0.,0.); 
00372   HepPoint3D hcenter(x_hc,y_hc,0.0);
00373 
00374   double m_c_perp(hcenter.perp());
00375   Hep3Vector m_c_unit((HepPoint3D)hcenter.unit());
00376   HepPoint3D IO = (-1) * charge * m_c_unit;
00377 
00378   dphi0  = fmod(IO.phi()+4*pi, 2*pi) - phi;
00379   IO_phi = fmod(IO.phi()+4*pi, 2*pi);
00380 
00381   if(dphi0 > pi) dphi0 -= 2*pi;
00382   else if(dphi0 < -pi) dphi0 += 2*pi; 
00383 
00384   phi_io[0] = -(1+charge)*pi/2 - charge*dphi0;
00385   phi_io[1] = phi_io[0]+1.5*pi;
00386 
00387 
00388   bool outFlag = false;
00389   //retrieve Mdc geometry information
00390   rCize1 = 0.1 * m_mdcGeomSvc->Layer(lay)->RCSiz1();  //mm -> cm 
00391   rCize2 = 0.1 * m_mdcGeomSvc->Layer(lay)->RCSiz2();  //mm -> cm
00392   rLay   = 0.1 * m_mdcGeomSvc->Layer(lay)->Radius();  //mm -> cm
00393   length = 0.1 * m_mdcGeomSvc->Layer(lay)->Length();  //mm -> cm
00394   //double halfLength=length/2.;
00395   //(conversion of the units mm(mc)->cm(rec))
00396   nCell  = m_mdcGeomSvc->Layer(lay)->NCell();
00397   phioffset = m_mdcGeomSvc->Layer(lay)->Offset();
00398   slant  = m_mdcGeomSvc->Layer(lay)->Slant();
00399   shift  = m_mdcGeomSvc->Layer(lay)->Shift();
00400   type   = m_mdcGeomSvc->Layer(lay)->Sup()->Type();
00401 
00402   m_crio[0] = rLay - rCize1 ; //radius of inner field wir ( beam pipe)
00403   m_crio[1] = rLay + rCize2 ; //radius of outer field wir ( beam pipe)
00404 
00405   int sign = -1; //assumed the first half circle
00406 
00407   Hep3Vector iocand[2];
00408   Hep3Vector cell_IO[2];
00409 
00410   for(int ii =0; ii<2; ii++){
00411     // By law of cosines to calculate the alpha(up and down field wir_Ge)
00412     double cos_alpha = (m_c_perp*m_c_perp + m_crio[ii]*m_crio[ii] - rho*rho)
00413       / ( 2 * m_c_perp * m_crio[ii] ); 
00414     if(fabs(cos_alpha)>1&&(ii==0)){
00415       cos_alpha = -1;
00416       outFlag=true;
00417     }
00418     if(fabs(cos_alpha)>1&&(ii==1)){
00419       cos_alpha = (m_c_perp*m_c_perp + m_crio[0]*m_crio[0] - rho*rho)
00420         / ( 2 * m_c_perp * m_crio[0] );
00421       C_alpha   = 2*pi - acos( cos_alpha);
00422     }else {
00423       C_alpha   = acos( cos_alpha );
00424     }
00425 
00426     iocand[ii] = m_c_unit;
00427     iocand[ii].rotateZ( charge*sign*C_alpha );
00428     iocand[ii] *= m_crio[ii];
00429 
00430     xx = iocand[ii].x() - x_hc ;
00431     yy = iocand[ii].y() - y_hc ;
00432 
00433     dphi[ii] = atan2(yy,xx) - phi0 - 0.5*pi*(1-charge);
00434     dphi[ii] = fmod( dphi[ii] + 8.0*pi,2*pi);
00435 
00436 
00437     if( dphi[ii] < phi_io[0] ) {
00438       dphi[ii] += 2*pi;
00439     }else if( dphi[ii] > phi_io[1] ){  //very very nausea
00440       dphi[ii] -= 2*pi;
00441     }
00442 
00443     cell_IO[ii] = Hel(piv, dr0, phi, ALPHA_loc, kappa,dz0, dphi[ii], tanl); 
00444 
00445     //cout<<" cell_IO["<<ii<<"] : "<<cell_IO[ii]<<endl;
00446     if( (cell_IO[ii].x()==0 ) && (cell_IO[ii].y()==0) && (cell_IO[ii].z()==0)) continue;
00447   }
00448   //if(((fabs(cell_IO[0].z())*10-halfLength)>-7.) && ((fabs(cell_IO[1].z())*10-halfLength)>-7.))return true; //Out sensitive area
00449 
00450   cellId_in = cellId_out = -1 ;
00451   phi_in  = cell_IO[0].phi();
00452   phi_in = fmod ( phi_in + 4 * pi, 2 * pi );
00453   phi_out = cell_IO[1].phi();
00454   phi_out = fmod ( phi_out + 4 * pi, 2 * pi );
00455   phi_bin = 2.0 * pi / nCell ;
00456 
00457   //determine the in/out cell id
00458   stphi[0] = shift * phi_bin * (0.5 - cell_IO[0].z()/length);
00459   stphi[1] = shift * phi_bin * (0.5 - cell_IO[1].z()/length);
00460   //stphi[0],stphi[1] to position fo z axis ,the angle of deflxsion in x_Y
00461 
00462   phioff[0] = phioffset + stphi[0];
00463   phioff[1] = phioffset + stphi[1];
00464 
00465   for(int kk = 0; kk<nCell ; kk++){
00466     //for stereo layer
00467     inlow  = phioff[0] + phi_bin*kk - phi_bin*0.5;
00468     inlow  = fmod( inlow + 4.0 * pi , 2.0 * pi);
00469     inup   = phioff[0] + phi_bin*kk + phi_bin*0.5;
00470     inup   = fmod( inup  + 4.0 * pi , 2.0 * pi);
00471     outlow = phioff[1] + phi_bin*kk - phi_bin*0.5;
00472     outlow = fmod( outlow + 4.0 * pi ,2.0 * pi);
00473     outup  = phioff[1] + phi_bin*kk + phi_bin*0.5;
00474     outup  = fmod( outup + 4.0 * pi , 2.0 * pi);
00475 
00476     if((phi_in>=inlow && phi_in<=inup))   cellId_in = kk;
00477     if((phi_out>=outlow&&phi_out<outup))  cellId_out = kk;
00478     if(inlow > inup ){
00479       if((phi_in>=inlow&&phi_in<2.0*pi)||(phi_in>=0.0&&phi_in<inup)) cellId_in = kk;
00480     }
00481     if(outlow>outup){
00482       if((phi_out>=outlow&&phi_out<=2.0*pi)||(phi_out>=0.0&&phi_out<outup)) cellId_out = kk;
00483     }
00484   }//end of nCell loop  
00485 
00486   return (cellId_in==cellId_out);
00487 }

double MdcCheckUtil::doca (  int  layer,
int  cell,
const MdcSWire *  sWire,
const HepVector  helix,
const HepSymMatrix  errMat
) 
 

double MdcCheckUtil::doca (  int  layer,
int  cell,
HepPoint3D  eastP,
HepPoint3D  westP,
const HepVector  helix,
const HepSymMatrix  errMat
) 
 

double MdcCheckUtil::doca (  int  layer,
int  cell,
const HepVector  helix,
const HepSymMatrix  errMat
) 
 

double MdcCheckUtil::doca (  int  layer,
int  cell,
const MdcSWire *  sWire,
const HepVector  helix,
const HepSymMatrix  errMat
) 
  

00232                                                                                                                           {
00233 
00234   HepVector helixPat = besPar2PatPar(helixBes);
00235   HepSymMatrix errMatPat= besErr2PatErr(errMatBes);
00236 
00237   return docaPatPar(layer, wire, sWire, helixPat, errMatPat);
00238   
00239 }

double MdcCheckUtil::doca (  int  layer,
int  cell,
HepPoint3D  eastP,
HepPoint3D  westP,
const HepVector  helix,
const HepSymMatrix  errMat
) 
  

00216                                                                                                                                         {
00217   // get doca by user referenced point and layer cell
00218   double sag = m_gm->Wire(layer,cell)->getSag();
00219   int id = m_gm->Wire(layer,cell)->Id();
00220   //const MdcSWire* m_mdcSWire  = new MdcSWire(m_mdcGeomSvc, id);
00221   const MdcSWire* m_mdcSWire  = new MdcSWire(eastP, westP, sag, id , cell);
00222 
00223   double mdoca = doca(layer, cell, m_mdcSWire, helixBes, errMatBes);
00224 
00225   delete m_mdcSWire;
00226   return mdoca;
00227 }

double MdcCheckUtil::doca (  int  layer,
int  cell,
const HepVector  helix,
const HepSymMatrix  errMat
) 
  

00201                                                                                                     {
00202   const MdcSWire* m_mdcSWire  = m_gm->Wire(layer,cell);
00203   return doca(layer, cell, m_mdcSWire, helixBes, errMatBes);
00204 }

double MdcCheckUtil::docaPatPar (  int  layer,
int  cell,
const MdcSWire *  sWire,
const HepVector  helixPat,
const HepSymMatrix  errMatPat
) 
 

double MdcCheckUtil::docaPatPar (  int  layer,
int  cell,
const HepVector  helixPat,
const HepSymMatrix  errMatPat
) 
 

double MdcCheckUtil::docaPatPar (  int  layer,
int  cell,
const MdcSWire *  sWire,
const HepVector  helixPat,
const HepSymMatrix  errMatPat
) 
  

00242                                                                                                                           {
00243 
00244   //-----cell track passed
00245   int cellId_in  = -1;
00246   int cellId_out = -1;
00247 
00248   cellTrackPassed(helix,layer,cellId_in,cellId_out);
00249   // std::cout<<" cellId in  "<<cellId_in<<" out "<<cellId_out
00250   //  <<" pass 1 cell "<<passedOneCell<<std::endl;
00251   //cout << "wire = " << wire << ", cellId_in = " << cellId_in << ", cellId_out = " << cellId_out << endl;
00252   if (m_passCellRequired &&(wire < cellId_in && wire > cellId_out)) return -999.;
00253 
00254   //-----helix trajectory
00255   HelixTraj* m_helixTraj = new HelixTraj(helix,errMat);
00256   const Trajectory* trajHelix    = dynamic_cast<const Trajectory*> (m_helixTraj);
00257 
00258   //-----pointOnHelix
00259   int innerOrOuter = 1;
00260   Hep3Vector cell_IO = pointOnHelixPatPar(helix,layer,innerOrOuter);
00261   double fltTrack = 0.1 * m_mdcGeomSvc -> Layer(layer)->Radius();
00262 
00263   //------wire trajectory
00264   const MdcSagTraj* m_wireTraj = sWire->getTraj();
00265   const Trajectory* trajWire  = dynamic_cast<const Trajectory*> (m_wireTraj);
00266   const HepPoint3D* start_In  = sWire->getEastPoint();
00267   //const HepPoint3D* stop_In   = sWire->getWestPoint();
00268   //std::cout<< " east -------- "<< start_In->x()<<","<<start_In->y()<<","<<start_In->z() <<std::endl;
00269   //std::cout<< " west -------- "<< stop_In->x()<<","<<stop_In->y()<<","<<stop_In->z() <<std::endl;
00270 
00271   //------calc poca
00272   double doca  = -999.;
00273   TrkPoca* m_trkPoca;
00274   double zWire = cell_IO.z();
00275   HepPoint3D pos_in(zWire,sWire->xWireDC(zWire),sWire->yWireDC(zWire)) ;
00276 
00277   double fltWire = sqrt( (pos_in.x()-start_In->x())*(pos_in.x()-start_In->x()) +
00278       (pos_in.y()-start_In->y())*(pos_in.y()-start_In->y()) +
00279       (pos_in.z()-start_In->z())*(pos_in.z()-start_In->z())  );
00280   m_trkPoca  = new TrkPoca(*trajHelix, fltTrack, *trajWire, fltWire);
00281   doca = m_trkPoca->doca();
00282 
00283   delete m_helixTraj;
00284   delete m_trkPoca;
00285 
00286   return doca;
00287 
00288 } //----------end of calculatng the doca ---------------------------------//

double MdcCheckUtil::docaPatPar (  int  layer,
int  cell,
const HepVector  helixPat,
const HepSymMatrix  errMatPat
) 
  

00208                                                                                                           {
00209   const MdcSWire* m_mdcSWire  = m_gm->Wire(layer,cell);
00210   return docaPatPar(layer, cell, m_mdcSWire, helixPat, errMatPat);
00211 }

void MdcCheckUtil::dumpRecMdcTrack (   ) 
 

void MdcCheckUtil::dumpRecMdcTrack (   ) 
  

00648                                   {
00649   IDataProviderSvc* eventSvc = NULL;
00650   Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
00651   if (NULL == eventSvc){
00652     std::cout<<__FILE__<<" Could not load EventDataSvc service"<<std::endl;
00653     return;
00654   }
00655 
00656   SmartDataPtr<RecMdcTrackCol> trackList(eventSvc,EventModel::Recon::RecMdcTrackCol);
00657   if (!trackList) return;
00658   std::cout<< "tksize = "<<trackList->size() << std::endl;//yzhang debug
00659   RecMdcTrackCol::iterator it = trackList->begin();
00660   for (;it!= trackList->end();it++){
00661     RecMdcTrack *tk = *it;
00662     std::cout<< "//====RecMdcTrack "<<tk->trackId()<<"====:" << std::endl;
00663     cout <<" d0 "<<tk->helix(0)
00664       <<" phi0 "<<tk->helix(1)
00665       <<" cpa "<<tk->helix(2)
00666       <<" z0 "<<tk->helix(3)
00667       <<" tanl "<<tk->helix(4)
00668       <<endl;
00669     std::cout<<" q "<<tk->charge() 
00670       <<" theta "<<tk->theta()
00671       <<" phi "<<tk->phi()
00672       <<" x0 "<<tk->x()
00673       <<" y0 "<<tk->y()
00674       <<" z0 "<<tk->z()
00675       <<" r0 "<<tk->r()
00676       <<endl;
00677     std::cout <<" p "<<tk->p()
00678       <<" pt "<<tk->pxy()
00679       <<" px "<<tk->px()
00680       <<" py "<<tk->py()
00681       <<" pz "<<tk->pz()
00682       <<endl;
00683     std::cout<<" tkStat "<<tk->stat()
00684       <<" chi2 "<<tk->chi2()
00685       <<" ndof "<<tk->ndof()
00686       <<" nhit "<<tk->getNhits()
00687       <<" nst "<<tk->nster()
00688       <<endl;
00689     //std::cout<< "errmat   " << std::endl;
00690     //for (int i=0; i<15; i++){ std::cout<< " "<<tk->err(i); }
00691     //std::cout<< "   " << std::endl;
00692 
00693     int nhits = tk->getVecHits().size();
00694     std::cout<<nhits <<" Hits: " << std::endl;
00695     for(int ii=0; ii <nhits ; ii++){
00696       Identifier id(tk->getVecHits()[ii]->getMdcId());
00697       int layer = MdcID::layer(id);
00698       int wire = MdcID::wire(id);
00699       cout<<"("<< layer <<","<<wire<<","<<tk->getVecHits()[ii]->getStat()
00700         <<",lr:"<<tk->getVecHits()[ii]->getFlagLR()<<") ";
00701     }//end of hit list
00702     std::cout << "  "<< std::endl;
00703   }//end of tk list
00704   std::cout << "  "<< std::endl;
00705 }

HepPoint3D MdcCheckUtil::Hel (  HepPoint3D  piv,
double  dr,
double  phi0,
double  Alpha_L,
double  kappa,
double  dz,
double  dphi,
double  tanl
) 
 

HepPoint3D MdcCheckUtil::Hel (  HepPoint3D  piv,
double  dr,
double  phi0,
double  Alpha_L,
double  kappa,
double  dz,
double  dphi,
double  tanl
) 
  

00489                                                                                                                                {
00490   double x = piv.x() + dr*cos(phi0) + (Alpha_L/kappa) * (cos(phi0) - cos(phi0+dphi));
00491   double y = piv.y() + dr*sin(phi0) + (Alpha_L/kappa) * (sin(phi0) - sin(phi0+dphi));
00492   double z = piv.z() + dz - (Alpha_L/kappa) * dphi * tanl;
00493   //cout<<"HepPoint3D(x, y, z) = "<<HepPoint3D(x, y, z)<<endl;
00494   if((x>-1000. && x<1000.) || (y >-1000. && y <1000. ) ||(z>-1000. && z<1000.)){
00495     return HepPoint3D(x, y, z);
00496   }else{
00497     return HepPoint3D(0,0,0);
00498   }
00499 }

Hep3Vector MdcCheckUtil::momentum (  const RecMdcTrack *  trk  ) 
 

Hep3Vector MdcCheckUtil::momentum (  const RecMdcTrack *  trk  ) 
  

00579                                                         {
00580   // double dr   = trk->helix(0);
00581   double fi0  = trk->helix(1);
00582   double cpa  = trk->helix(2);
00583   double tanl = trk->helix(4);
00584 
00585   double pxy = 0.;
00586   if(cpa != 0) pxy = 1/fabs(cpa);
00587 
00588   double px = pxy * (-sin(fi0));
00589   double py = pxy * cos(fi0);
00590   double pz = pxy * tanl;
00591 
00592   Hep3Vector p;
00593   p.set(px, py, pz); // MeV/c
00594   return p;
00595 }

int MdcCheckUtil::nLayerTrackPassed (  const double  helix[5]  ) 
 

int MdcCheckUtil::nLayerTrackPassed (  const HepVector  helix  ) 
 

int MdcCheckUtil::nLayerTrackPassed (  const double  helix[5]  ) 
  

00067                                                            {
00068 
00069   HepVector helixBesParam(5,0);
00070   for(int i=0; i<5; ++i) helixBesParam[i] = helixBes[i];
00071 
00072   return nLayerTrackPassed(helixBesParam);
00073 }

int MdcCheckUtil::nLayerTrackPassed (  const HepVector  helix  ) 
  

00076                                                            {
00077   HepVector helix= besPar2PatPar(helixBes);
00078 
00079   int nLayer = 0;
00080 
00081   for(unsigned iLayer=0; iLayer<43; iLayer++){
00082     //flightLength is the path length of track in the x-y plane
00083     //guess flightLength by the radius in middle of the wire.
00084     double rMidLayer = m_mdcGeomSvc->Layer(iLayer)->Radius();
00085     double flightLength = rMidLayer;
00086 
00087     HepPoint3D pivot(0.,0.,0.);
00088     double dz = helix[3];
00089     double c = CLHEP::c_light * 100.; //unit from m/s to cm/s
00090     double alpha = 1/(c * Bz);//~333.567
00091     double kappa = helix[2];
00092     double rc = (-1.)*alpha/kappa;
00093     //std::cout<<"MdcCheckUtil alpha   "<<alpha<<std::endl;
00094     double tanl = helix[4];
00095     double phi0 = helix[1];
00096     double phi = flightLength/rc + phi0;//turning angle
00097     double z = pivot.z() + dz - (alpha/kappa) * tanl * phi;
00098 
00099     double layerHalfLength = m_mdcGeomSvc->Layer(iLayer)->Length()/2.;
00100 
00101     //std::cout<<"MdcCheckUtil length  "<<layerHalfLength<<std::endl;
00102 
00103     if (fabs(z) < fabs(layerHalfLength)) ++nLayer;
00104   }
00105 
00106   return nLayer;
00107 }

double MdcCheckUtil::p_cms (  HepVector  helix,
int  runNo,
double  mass
) 
 

double MdcCheckUtil::p_cms (  HepVector  helix,
int  runNo,
double  mass
) 
  

00547                                                                  {
00548   HepLorentzVector p4;
00549   p4.setPx(- sin(helix[1]) / fabs(helix[2]));
00550   p4.setPy(cos(helix[1]) / fabs(helix[2]));
00551   p4.setPz(helix[4] / fabs(helix[2]));
00552   double p3 = p4.mag();
00553   mass = 0.000511;
00554   p4.setE(sqrt(p3 * p3 + mass * mass));
00555 
00556   double ecm;
00557   if (runNo > 9815) {
00558     ecm = 3.097;
00559   }else{
00560     ecm = 3.68632;
00561   }
00562   double zboost = 0.0075;
00563   //HepLorentzVector psip(0.011 * 3.68632, 0, 0.0075, 3.68632);
00564   HepLorentzVector psip(0.011 * ecm, 0, zboost, ecm);
00565   //cout << setw(15) << ecm << setw(15) << zboost << endl;
00566   Hep3Vector boostv = psip.boostVector();
00567 
00568   //std::cout<<__FILE__<<" boostv "<<boostv<<  std::endl;
00569   p4.boost(- boostv);
00570  
00571   //std::cout<<__FILE__<<" p4 "<<p4<<  std::endl;
00572   double p_cms = p4.rho();
00573   //phicms = p4.phi();
00574   //p4.boost(boostv);
00575 
00576   return p_cms;
00577 }

HepSymMatrix MdcCheckUtil::patRecErr2BesErr (  const HepSymMatrix &  err  ) 
 

HepSymMatrix MdcCheckUtil::patRecErr2BesErr (  const HepSymMatrix &  err  ) 
  

00127                                                                   {
00128   //error matrix
00129   //std::cout<<" err1  "<<err<<" "<<err.num_row()<<std::endl;
00130   //V(Y)=S * V(X) * ST , mS = S , mVy = V(Y) , err() = V(X)
00131   //int n = err.num_row();
00132   HepSymMatrix mS(err.num_row(),0);
00133   mS[0][0]=-1.;//dr0=-d0
00134   mS[1][1]=1.;
00135   mS[2][2]=Bz/-333.567;//pxy -> cpa
00136   mS[3][3]=1.;
00137   mS[4][4]=1.;
00138   HepSymMatrix mVy= err.similarity(mS);
00139   //std::cout<<" err2  "<<n<<" "<<mVy<<std::endl;
00140   return mVy;
00141 }

HepVector MdcCheckUtil::patRecPar2BesPar (  const HepVector &  helixPar  ) 
 

HepVector MdcCheckUtil::patRecPar2BesPar (  const HepVector &  helixPar  ) 
  

00110                                                                  {
00111   HepVector helix(5,0);
00112   double d0 = -helixPar[0];    //cm
00113   double phi0 = helixPar[1]+ CLHEP::halfpi;
00114   double omega = Bz*helixPar[2]/-333.567;
00115   double z0 = helixPar[3];    //cm
00116   double tanl = helixPar[4];
00117   helix[0] = d0;
00118   helix[1] = phi0;
00119   helix[2] = omega;
00120   helix[3] = z0;
00121   helix[4] = tanl;
00122   //std::cout<<"helix   "<<helix<<std::endl;
00123   return helix;
00124 }

HepPoint3D MdcCheckUtil::pointOnHelix (  const HepVector  helixPar,
int  lay,
int  innerOrOuter
) 
 

HepPoint3D MdcCheckUtil::pointOnHelix (  const HepVector  helixPar,
int  lay,
int  innerOrOuter
) 
  

00145                                                                                           {
00146   HepVector helix= besPar2PatPar(helixBes);
00147   return pointOnHelixPatPar(helix, layer, innerOrOuter);
00148 }

HepPoint3D MdcCheckUtil::pointOnHelixPatPar (  const HepVector  helixPat,
int  lay,
int  innerOrOuter
) 
 

HepPoint3D MdcCheckUtil::pointOnHelixPatPar (  const HepVector  helixPat,
int  lay,
int  innerOrOuter
) 
  

00152                                                                                              {
00153   
00154   double rInner,rOuter;
00155   //retrieve Mdc geometry information
00156   double rCize1 =0.1* m_mdcGeomSvc->Layer(layer)->RCSiz1();  //mm -> cm 
00157   double rCize2 =0.1* m_mdcGeomSvc->Layer(layer)->RCSiz2();  //mm -> cm
00158   double rLay   =0.1* m_mdcGeomSvc->Layer(layer)->Radius();  //mm -> cm
00159 
00160   //(conversion of the units mm(mc)->cm(rec))
00161   rInner = rLay - rCize1 ; //radius of inner field wire
00162   rOuter = rLay + rCize2 ; //radius of outer field wire
00163 
00164   //int sign = -1; //assumed the first half circle
00165   HepPoint3D piv(0.,0.,0.);
00166   double r;
00167   if (innerOrOuter) r = rInner;
00168   else r = rOuter;
00169 
00170   double d0 = helix[0];
00171   double phi0 = helix[1];
00172   double omega = helix[2];
00173   double z0 = helix[3];
00174   double tanl = helix[4];
00175 
00176   double rc; 
00177   if (abs(omega) <Constants::epsilon) rc = 9999.;
00178   else rc = 1./omega;
00179   double xc     = piv.x() + ( d0 + rc) * cos(phi0);
00180   double yc     = piv.y() + ( d0 + rc ) * sin(phi0);
00181   HepPoint3D helixCenter(xc,yc,0.0);
00182   rc = sqrt(xc*xc + yc*yc);
00183   double a,b,c;
00184   a = r;
00185   b = d0 + rc;
00186   c = rc;
00187   double dphi = acos((a*a-b*b-c*c)/(-2.*b*c)); 
00188   double fltlen = dphi * rc;
00189   double phi = phi0 * omega * fltlen;
00190   double x = piv.x()+d0*sin(phi) - (rc+d0)*sin(phi0);
00191   double y = piv.y()+d0*cos(phi) + (rc+d0)*cos(phi0);
00192   double z = piv.z()+ z0 + fltlen*tanl;
00193   //std::cout<<" in Hel  "<<helix<<" "<<r<<" "<<rc<<" "<<dphi<<std::endl;
00194   //cout<<"HepPoint3D(x, y, z) = "<<HepPoint3D(x, y, z)<<endl;
00195   return HepPoint3D(x, y, z);
00196 }

double MdcCheckUtil::probab (  const int &  ndof,
const double &  chisq
) 
 

double MdcCheckUtil::probab (  const int &  ndof,
const double &  chisq
) 
  

00597                                                                {
00598 
00599   //constants
00600   double srtopi=2.0/sqrt(2.0*M_PI);
00601   double upl=100.0;
00602 
00603   double prob=0.0;
00604   if(ndof<=0) {return prob;}
00605   if(chisq<0.0) {return prob;}
00606   if(ndof<=60) {
00607     //full treatment
00608     if(chisq>upl) {return prob;}
00609     double sum=exp(-0.5*chisq);
00610     double term=sum;
00611 
00612     int m=ndof/2;
00613     if(2*m==ndof){
00614       if(m==1){return sum;}
00615       for(int i=2; i<=m;i++){
00616         term=0.5*term*chisq/(i-1);
00617         sum+=term;
00618       }//(int i=2; i<=m)
00619       return sum;
00620       //even
00621 
00622     }else{
00623       //odd
00624       double srty=sqrt(chisq);
00625       double temp=srty/M_SQRT2;
00626       prob=erfc(temp);
00627       if(ndof==1) {return prob;}
00628       if(ndof==3) {return (srtopi*srty*sum+prob);}
00629       m=m-1;
00630       for(int i=1; i<=m; i++){
00631         term=term*chisq/(2*i+1);
00632         sum+=term;
00633       }//(int i=1; i<=m; i++)
00634       return (srtopi*srty*sum+prob);
00635 
00636     }//(2*m==ndof)
00637 
00638   }else{
00639     //asymtotic Gaussian approx
00640     double srty=sqrt(chisq)-sqrt(ndof-0.5);
00641     if(srty<12.0) {prob=0.5*erfc(srty);};
00642 
00643   }//ndof<30
00644 
00645   return prob;
00646 }//endof probab

void MdcCheckUtil::setPassCellRequired (  bool  pass  )  [inline]
  

00045 {m_passCellRequired = pass;}

void MdcCheckUtil::setPassCellRequired (  bool  pass  )  [inline]
  

00045 {m_passCellRequired = pass;}

Member Data Documentation

double MdcCheckUtil::Bz [private]
 

int MdcCheckUtil::m_debug [private]
 

const MdcDetector* MdcCheckUtil::m_gm [private]
 

const MdcDetector* MdcCheckUtil::m_gm [private]
 

MdcGeomSvc* MdcCheckUtil::m_mdcGeomSvc [private]
 

MdcGeomSvc* MdcCheckUtil::m_mdcGeomSvc [private]
 

bool MdcCheckUtil::m_passCellRequired [private]
 

IMagneticFieldSvc* MdcCheckUtil::m_pIMF [private]
 

IMagneticFieldSvc* MdcCheckUtil::m_pIMF [private]
 

The documentation for this class was generated from the following files:
Generated on Wed Feb 2 16:24:17 2011 for BOSS6.5.5 by  doxygen 1.3.9.1