MdcUtilitySvc Class Reference

#include <MdcUtilitySvc.h>

Inheritance diagram for MdcUtilitySvc:

Public Member Functions
	MdcUtilitySvc (const std::string &name, ISvcLocator *svcloc)
	~MdcUtilitySvc ()
virtual StatusCode	initialize ()
virtual StatusCode	finalize ()
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
MdcDigiVec	getMdcDigiVec () const
void	getHelixOfMcParticle (const Event::McParticle *mcParticle, Helix &helix)
float	getChargeOfMcParticle (const Event::McParticle *mcParticle)
void	getMomPosOfMcParticle (const Event::McParticle *mcParticle, HepVector3D &pos, HepVector3D &mom)
void	getMdcDigiOnMcParticle (int trackIndex, const std::map< int, std::map< MdcDigi *, Event::MdcMcHit * > > mdcMCAssociation, MdcDigiVec &mdcDigiInput, MdcDigiVec &mdcDigiAssociated)
void	getMdcMCAssoiciation (int trackIndex, const std::vector< MdcDigi * > mdcDigiVecInput, std::map< int, std::map< MdcDigi *, Event::MdcMcHit * > > &mdcMCAssociation)
	Get association of MdcDigi and MdcMcHit according to track id.

Detailed Description

Definition at line 40 of file MdcUtilitySvc.h.

Constructor & Destructor Documentation

MdcUtilitySvc()

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

Definition at line 24 of file MdcUtilitySvc.cxx.

                                                                     :
  base_class (name, svcloc) {
    declareProperty("debug",        m_debug = 0);
  }

~MdcUtilitySvc()

MdcUtilitySvc::~MdcUtilitySvc

(

)

Definition at line 29 of file MdcUtilitySvc.cxx.

                             {
}

Member Function Documentation

besErr2PatErr()

HepSymMatrix MdcUtilitySvc::besErr2PatErr

(

const HepSymMatrix &

err

)

const

Definition at line 777 of file MdcUtilitySvc.cxx.

                                                                      {
  //error matrix
  //std::cout<<" err1  "<<err<<" "<<err.num_row()<<std::endl;
  //V(Y)=S * V(X) * ST , mS = S , mVy = V(Y) , err() = V(X)
  //int n = err.num_row();
  HepSymMatrix mS(err.num_row(),0);
  mS[0][0]=-1.;//dr0=-d0
  mS[1][1]=1.;
  mS[2][2]=Bz()/-333.567;//pxy -> cpa
  mS[3][3]=1.;
  mS[4][4]=1.;
  HepSymMatrix mVy= err.similarity(mS);
  //std::cout<<" err2  "<<n<<" "<<mVy<<std::endl;
  return mVy;
}

Referenced by doca(), doca(), and doca().

besPar2PatPar()

HepVector MdcUtilitySvc::besPar2PatPar

(

const HepVector &

helixPar

)

const

Definition at line 761 of file MdcUtilitySvc.cxx.

                                                                     {
  HepVector helix(5,0);
  double d0 = -helixPar[0];    //cm
  double phi0 = helixPar[1]+ CLHEP::halfpi;
  double omega = Bz()*helixPar[2]/-333.567;
  double z0 = helixPar[3];    //cm
  double tanl = helixPar[4];
  helix[0] = d0;
  helix[1] = phi0;
  helix[2] = omega;
  helix[3] = z0;
  helix[4] = tanl;
  return helix;
}

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

cellTrackPassed()

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

Definition at line 599 of file MdcUtilitySvc.cxx.

                                                                                                   {
  int charge,type,nCell;
  double dr0,phi0,kappa,dz0,tanl;
  double ALPHA_loc,rho,phi,cosphi0,sinphi0,x_hc,y_hc,z_hc;
  double dphi0,IO_phi,C_alpha,xx,yy;
  double inlow,inup,outlow,outup,phi_in,phi_out,phi_bin;
  double rCize1,rCize2,rLay,length,phioffset,slant,shift;
  double m_crio[2], phi_io[2], stphi[2],phioff[2],dphi[2];
 
  dr0   = helixBes[0];
  phi0  = helixBes[1];
  kappa = helixBes[2];
  dz0   = helixBes[3];
  tanl  = helixBes[4];
 
  ALPHA_loc = 1000/(2.99792458*Bz()); //magnetic field const
  charge    = ( kappa >=0 ) ? 1 : -1 ;
  rho       = ALPHA_loc/kappa ;
  double pi = Constants::pi;
  phi      = fmod(phi0 + 4*pi , 2*pi);
  cosphi0  = cos(phi);
  sinphi0  = (1.0 - cosphi0 ) * (1.0 + cosphi0 );
  sinphi0  = sqrt(( sinphi0 > 0.) ? sinphi0 : 0.);
  if( phi > pi ) sinphi0 = -sinphi0 ;
 
  x_hc     = 0. + ( dr0 + rho ) * cosphi0;
  y_hc     = 0. + ( dr0 + rho ) * sinphi0;
  z_hc     = 0. + dz0;
 
 
  HepPoint3D piv(0.,0.,0.);
  HepPoint3D hcenter(x_hc,y_hc,0.0);
 
  double m_c_perp(hcenter.perp());
  Hep3Vector m_c_unit((HepPoint3D)hcenter.unit());
  HepPoint3D IO = (-1) * charge * m_c_unit;
 
  dphi0  = fmod(IO.phi()+4*pi, 2*pi) - phi;
  IO_phi = fmod(IO.phi()+4*pi, 2*pi);
 
  if(dphi0 > pi) dphi0 -= 2*pi;
  else if(dphi0 < -pi) dphi0 += 2*pi;
 
  phi_io[0] = -(1+charge)*pi/2 - charge*dphi0;
  phi_io[1] = phi_io[0]+1.5*pi;
 
 
  bool outFlag = false;
  //retrieve Mdc geometry information
  rCize1 = 0.1 * m_mdcGeomSvc->Layer(lay)->RCSiz1();  //mm -> cm
  rCize2 = 0.1 * m_mdcGeomSvc->Layer(lay)->RCSiz2();  //mm -> cm
  rLay   = 0.1 * m_mdcGeomSvc->Layer(lay)->Radius();  //mm -> cm
  length = 0.1 * m_mdcGeomSvc->Layer(lay)->Length();  //mm -> cm
  //double halfLength=length/2.;
  //(conversion of the units mm(mc)->cm(rec))
  nCell  = m_mdcGeomSvc->Layer(lay)->NCell();
  phioffset = m_mdcGeomSvc->Layer(lay)->Offset();
  slant  = m_mdcGeomSvc->Layer(lay)->Slant();
  shift  = m_mdcGeomSvc->Layer(lay)->Shift();
  type   = m_mdcGeomSvc->Layer(lay)->Sup()->Type();
 
  m_crio[0] = rLay - rCize1 ; //radius of inner field wir ( beam pipe)
  m_crio[1] = rLay + rCize2 ; //radius of outer field wir ( beam pipe)
 
  int sign = -1; //assumed the first half circle
 
  Hep3Vector iocand[2];
  Hep3Vector cell_IO[2];
 
  for(int ii =0; ii<2; ii++){
    // By law of cosines to calculate the alpha(up and down field wir_Ge)
    double cos_alpha = (m_c_perp*m_c_perp + m_crio[ii]*m_crio[ii] - rho*rho)
      / ( 2 * m_c_perp * m_crio[ii] );
    if(fabs(cos_alpha)>1&&(ii==0)){
      cos_alpha = -1;
      outFlag=true;
    }
    if(fabs(cos_alpha)>1&&(ii==1)){
      cos_alpha = (m_c_perp*m_c_perp + m_crio[0]*m_crio[0] - rho*rho)
  / ( 2 * m_c_perp * m_crio[0] );
      C_alpha   = 2*pi - acos( cos_alpha);
    }else {
      C_alpha   = acos( cos_alpha );
    }
 
    iocand[ii] = m_c_unit;
    iocand[ii].rotateZ( charge*sign*C_alpha );
    iocand[ii] *= m_crio[ii];
 
    xx = iocand[ii].x() - x_hc ;
    yy = iocand[ii].y() - y_hc ;
 
    dphi[ii] = atan2(yy,xx) - phi0 - 0.5*pi*(1-charge);
    dphi[ii] = fmod( dphi[ii] + 8.0*pi,2*pi);
 
 
    if( dphi[ii] < phi_io[0] ) {
      dphi[ii] += 2*pi;
    }else if( dphi[ii] > phi_io[1] ){  //very very nausea
      dphi[ii] -= 2*pi;
    }
 
    cell_IO[ii] = Hel(piv, dr0, phi, ALPHA_loc, kappa,dz0, dphi[ii], tanl);
 
    //cout<<" cell_IO["<<ii<<"] : "<<cell_IO[ii]<<endl;
    if( (cell_IO[ii].x()==0 ) && (cell_IO[ii].y()==0) && (cell_IO[ii].z()==0)) continue;
  }
  //if(((fabs(cell_IO[0].z())*10-halfLength)>-7.) && ((fabs(cell_IO[1].z())*10-halfLength)>-7.))return true; //Out sensitive area
 
  cellId_in = cellId_out = -1 ;
  phi_in  = cell_IO[0].phi();
  phi_in = fmod ( phi_in + 4 * pi, 2 * pi );
  phi_out = cell_IO[1].phi();
  phi_out = fmod ( phi_out + 4 * pi, 2 * pi );
  phi_bin = 2.0 * pi / nCell ;
 
  //determine the in/out cell id
  stphi[0] = shift * phi_bin * (0.5 - cell_IO[0].z()/length);
  stphi[1] = shift * phi_bin * (0.5 - cell_IO[1].z()/length);
  //stphi[0],stphi[1] to position fo z axis ,the angle of deflxsion in x_Y
 
  phioff[0] = phioffset + stphi[0];
  phioff[1] = phioffset + stphi[1];
 
  for(int kk = 0; kk<nCell ; kk++){
    //for stereo layer
    inlow  = phioff[0] + phi_bin*kk - phi_bin*0.5;
    inlow  = fmod( inlow + 4.0 * pi , 2.0 * pi);
    inup   = phioff[0] + phi_bin*kk + phi_bin*0.5;
    inup   = fmod( inup  + 4.0 * pi , 2.0 * pi);
    outlow = phioff[1] + phi_bin*kk - phi_bin*0.5;
    outlow = fmod( outlow + 4.0 * pi ,2.0 * pi);
    outup  = phioff[1] + phi_bin*kk + phi_bin*0.5;
    outup  = fmod( outup + 4.0 * pi , 2.0 * pi);
 
    if((phi_in>=inlow && phi_in<=inup))   cellId_in = kk;
    if((phi_out>=outlow&&phi_out<outup))  cellId_out = kk;
    if(inlow > inup ){
      if((phi_in>=inlow&&phi_in<2.0*pi)||(phi_in>=0.0&&phi_in<inup)) cellId_in = kk;
    }
    if(outlow>outup){
      if((phi_out>=outlow&&phi_out<=2.0*pi)||(phi_out>=0.0&&phi_out<outup)) cellId_out = kk;
    }
  }//end of nCell loop
 
  return (cellId_in==cellId_out);
}

Referenced by cellTrackPassedPatPar().

cellTrackPassedByPhi()

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

Definition at line 552 of file MdcUtilitySvc.cxx.

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

cellTrackPassedByPhiPatPar()

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

Definition at line 561 of file MdcUtilitySvc.cxx.

                                                                                                                     {
  int nCell = m_mdcGeomSvc->Layer(layer)->NCell();
 
  double dPhiz = (m_mdcGeomSvc->Wire(layer,0)->Forward().phi() - m_mdcGeomSvc->Wire(layer,0)->Backward().phi())*0.5;
  double rEnd = m_mdcGeomSvc->Wire(layer,0)->Backward().rho()/10.;//mm2cm
  double rMid = rEnd * cos(dPhiz);
  //std::cout<<"( cell "<<0<<" westPphi "<<m_mdcGeomSvc->Wire(layer,0)->Forward().phi() <<" eastPphi "
  //<<m_mdcGeomSvc->Wire(layer,0)->Backward().phi()<<" twist  "<<dPhiz<<" rend "<<rEnd<<std::endl;
  double fltLenIn = rMid;
  double phiIn = helixPat[1] + helixPat[2]*fltLenIn;//phi0 + omega * L
 
  double phiEPOffset = m_mdcGeomSvc->Wire(layer,0)->Backward().phi();//east.phi()= BackWard.phi
  BesAngle tmp(phiIn - phiEPOffset);
  if(m_debug) std::cout<<"cellTrackPassed  nCell  "<<nCell<<" layer "<<layer<<" fltLenIn  "<<fltLenIn<<" phiEPOffset "<<phiEPOffset<<std::endl;
  //BesAngle tmp(phiIn - layPtr->phiEPOffset());
  int wlow = (int)floor(nCell * tmp.rad() / twopi );
  int wbig = (int)ceil(nCell * tmp.rad() / twopi );
  if (tmp == 0 ){ wlow = -1; wbig = 1; }
 
  if ((wlow%nCell)< 0) wlow += nCell;
  cellId_in  = wlow;
 
  if ((wbig%nCell)< 0) wbig += nCell;
  cellId_out = wbig;
 
  bool passedOneCell = (cellId_in == cellId_out);
 
  return passedOneCell;//pass more than one cell
}

Referenced by cellTrackPassedByPhi(), and docaPatPar().

cellTrackPassedPatPar()

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

Definition at line 592 of file MdcUtilitySvc.cxx.

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

doca() [1/3]

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

Definition at line 390 of file MdcUtilitySvc.cxx.

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

Referenced by docaPatPar(), docaPatPar(), docaPatPar(), and EsTimeAlg::execute().

doca() [2/3]

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

Definition at line 408 of file MdcUtilitySvc.cxx.

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

doca() [3/3]

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

Definition at line 399 of file MdcUtilitySvc.cxx.

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

docaPatPar() [1/3]

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

Definition at line 417 of file MdcUtilitySvc.cxx.

                                                                                                                                                 {
 
  const MdcGeoWire *geowir = m_mdcGeomSvc->Wire(layer,cell);
  int id = geowir->Id();
  double sag = 0.;
  if(doSag) sag = geowir->Sag()/10.;//mm2cm
  HepPoint3D eastP = geowir->Backward()/10.0;//mm2cm
  HepPoint3D westP = geowir->Forward() /10.0;//mm2cm
  const MdcSWire* sWire = new MdcSWire(eastP, westP, sag, id , cell);
 
  double doca = docaPatPar(layer, cell, sWire, helixPat, errMat, passCellRequired);//call 6.
 
  delete sWire;
  return doca;
}

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

docaPatPar() [2/3]

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

Definition at line 451 of file MdcUtilitySvc.cxx.

                                                                                                                                                           {
 
  if(m_debug) std::cout<<" helixPat  "<<helixPat<<std::endl;
  if(m_debug) std::cout<<" layer "<<layer<<" cell "<<cell<<std::endl;
  //-----cell track passed
  int cellId_in  = -1;
  int cellId_out = -1;
 
  //cellTrackPassedPatPar(helixPat,layer,cellId_in,cellId_out);//yzhang FIXME 2012-07-11
  cellTrackPassedByPhiPatPar(helixPat,layer,cellId_in,cellId_out);//yzhang FIXME 2012-07-11
 
  if(m_debug) {
    std::cout<<" cellId in  "<<cellId_in<<" out "<<cellId_out <<std::endl;
    cout << "cell = " << cell << ", cellId_in = " << cellId_in << ", cellId_out = " << cellId_out << endl;
  }
  if (passCellRequired &&(cell < cellId_in && cell > cellId_out)) return -999.;
 
  //-----helix trajectory
  HelixTraj* helixTraj = new HelixTraj(helixPat,errMat);
  const Trajectory* trajHelix    = dynamic_cast<const Trajectory*> (helixTraj);
 
  //-----pointOnHelix
  int innerOrOuter = 1;
  Hep3Vector cell_IO = pointOnHelixPatPar(helixPat,layer,innerOrOuter);
  double fltTrack = 0.1 * m_mdcGeomSvc -> Layer(layer)->Radius();
  if(m_debug) {
    std::cout<<" cell_IO "<<cell_IO<<std::endl;
    std::cout<<" fltTrack "<<fltTrack<<std::endl;
  }
 
  //------wire trajectory
  const MdcSagTraj* m_wireTraj = sWire->getTraj();
  const Trajectory* trajWire  = dynamic_cast<const Trajectory*> (m_wireTraj);
  const HepPoint3D* start_In  = sWire->getEastPoint();
  //const HepPoint3D* stop_In   = sWire->getWestPoint();
  if(m_debug){
    std::cout<<" sag "<<m_wireTraj->sag()<<std::endl;
    std::cout<< " east -------- "<< start_In->x()<<","<<start_In->y()<<","<<start_In->z()<<std::endl;
  }
  //std::cout<< " west -------- "<< stop_In->x()<<","<<stop_In->y()<<","<<stop_In->z() <<std::endl;
 
  //------calc poca
  double doca  = -999.;
  TrkPoca* trkPoca;
  double zWire = cell_IO.z();
  //HepPoint3D pos_in(zWire,sWire->xWireDC(zWire),sWire->yWireDC(zWire)) ;
  HepPoint3D pos_in(sWire->xWireDC(zWire),sWire->yWireDC(zWire),zWire) ;
  if(m_debug) std::cout<<" zWire  "<<zWire<<" zEndDC "<<sWire->zEndDC()<<std::endl;
  //if(m_debug) std::cout<<"pos_in "<<pos_in<<" fltWire  "<<fltWire<<std::endl;
 
  double fltWire = sqrt( (pos_in.x()-start_In->x())*(pos_in.x()-start_In->x()) +
      (pos_in.y()-start_In->y())*(pos_in.y()-start_In->y()) +
      (pos_in.z()-start_In->z())*(pos_in.z()-start_In->z())  );
  trkPoca  = new TrkPoca(*trajHelix, fltTrack, *trajWire, fltWire);
  doca = trkPoca->doca();
 
  double hitLen   = trkPoca->flt2();
  double startLen = trajWire->lowRange() - 5.;
  double endLen   = trajWire->hiRange()  + 5.;
  if(hitLen < startLen || hitLen > endLen) {
    if(m_debug) std::cout<<"WARNING MdcUtilitySvc::docaPatPar()  isBeyondEndflange! hitLen="<<hitLen <<" startLen="<<startLen<<" endLen "<<endLen<<std::endl;
    doca = -998.;
  }
  //std::cout<<" hitLen  "<<hitLen<<" startLen  "<<startLen<<" endLen  "<<endLen <<" doca "<<doca<<std::endl;
 
  if(m_debug) std::cout<<" doca  "<<doca<<std::endl;
  delete helixTraj;
  delete trkPoca;
 
  return doca;
 
} //----------end of calculatng the doca ---------------------------------//

docaPatPar() [3/3]

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

Definition at line 435 of file MdcUtilitySvc.cxx.

                                                                                                                                                                                  {
 
  const MdcGeoWire *geowir = m_mdcGeomSvc->Wire(layer,cell);
  int id = geowir->Id();
  double sag = 0.;
  if(doSag) sag = geowir->Sag()/10.;//mm2cm
  const MdcSWire* sWire = new MdcSWire(eastP, westP, sag, id , cell);//cm
 
  double doca = docaPatPar(layer, cell, sWire, helixPat, errMat, passCellRequired);//call 6.
 
  delete sWire;
  return doca;
}

finalize()

StatusCode MdcUtilitySvc::finalize ( )

virtual

Definition at line 70 of file MdcUtilitySvc.cxx.

                                  {
  MsgStream log(messageService(), name());
  log << MSG::INFO << "MdcUtilitySvc::finalize()" << endreq;
 
  return StatusCode::SUCCESS;
}

getChargeOfMcParticle()

float MdcUtilitySvc::getChargeOfMcParticle

(

const Event::McParticle *

mcParticle

)

Definition at line 895 of file MdcUtilitySvc.cxx.

                                                                           {
  // get PartPropSvc
  IPartPropSvc* partPropSvc;
  static const bool CREATEIFNOTTHERE(true);
  StatusCode sc=service("PartPropSvc",partPropSvc,CREATEIFNOTTHERE);
  if(!sc.isSuccess()||nullptr==partPropSvc){
    cout<< " Could not initialize Particle Properties Service" << endl;
  }
  HepPDT::ParticleDataTable* particleTable = partPropSvc->PDT();
  if(particleTable->particle(mcParticle->particleProperty())){
    return particleTable->particle(mcParticle->particleProperty())->charge();
  }
  return -1e6;
}

Referenced by getHelixOfMcParticle().

getHelixOfMcParticle()

void MdcUtilitySvc::getHelixOfMcParticle	(	const Event::McParticle *	mcParticle,
		Helix &	helix )

Definition at line 920 of file MdcUtilitySvc.cxx.

                                                                                      {
  if(nullptr==mcParticle) return;
  Hep3Vector pos,mom;
  pos.setX(mcParticle->initialPosition().x());
  pos.setY(mcParticle->initialPosition().y());
  pos.setZ(mcParticle->initialPosition().z());
  mom.setX(mcParticle->initialFourMomentum().px());
  mom.setY(mcParticle->initialFourMomentum().py());
  mom.setZ(mcParticle->initialFourMomentum().pz());
 
  Helix helixTemp(pos,mom,
      getChargeOfMcParticle(mcParticle));
  helixTemp.pivot(HepPoint3D(0,0,0));
  helix.set(helixTemp.pivot(),helixTemp.a(),helixTemp.Ea());
}

getMdcDigiOnMcParticle()

void MdcUtilitySvc::getMdcDigiOnMcParticle	(	int	trackIndex,
		const std::map< int, std::map< MdcDigi *, Event::MdcMcHit * > >	mdcMCAssociation,
		MdcDigiVec &	mdcDigiInput,
		MdcDigiVec &	mdcDigiAssociated )

Definition at line 936 of file MdcUtilitySvc.cxx.

                                                            {
  MdcDigiVec::iterator iter=mdcDigiInput.begin();
  //auto iter=mdcDigiVecInput.begin();
  for(;iter!=mdcDigiInput.end();iter++){
    std::map<int, std::map<MdcDigi*,MdcMcHit*> >::const_iterator itMap0=
      mdcMCAssociation.find(trackIndex);
    if(itMap0!=mdcMCAssociation.end()){
      std::map<MdcDigi*,MdcMcHit*>::const_iterator itMap=
        (*itMap0).second.find(*iter);
      MdcMcHit* mcHit;
      if(itMap!=(*itMap0).second.end()){ mcHit=itMap->second; }
      if( (trackIndex==mcHit->getTrackIndex())
          && (mcHit->getCreatorProcess()=="Generator"||
            mcHit->getCreatorProcess()=="Decay") ){
        mdcDigiAssociated.push_back(*iter);
      }
    }
  }
}

getMdcDigiVec()

vector< MdcDigi * > MdcUtilitySvc::getMdcDigiVec

(

)

const

Definition at line 77 of file MdcUtilitySvc.cxx.

{
  uint32_t getDigiFlag = 0;
 
  MdcDigiVec mdcDigiVec = m_RawDataProviderSvc->getMdcDigiVec(getDigiFlag);
  cout<<"MdcUtilitySvc::getMdcDigiVec() get "<<mdcDigiVec.size()<<" MDC digis from RawDataProviderSvc"<<endl;
 
  return mdcDigiVec;
}

getMdcMCAssoiciation()

void MdcUtilitySvc::getMdcMCAssoiciation	(	int	trackIndex,
		const std::vector< MdcDigi * >	mdcDigiVecInput,
		std::map< int, std::map< MdcDigi *, Event::MdcMcHit * > > &	mdcMCAssociation )

Get association of MdcDigi and MdcMcHit according to track id.

—get MDC MC hits and make a map vs wire id

loop over MdcDigi

Definition at line 958 of file MdcUtilitySvc.cxx.

                                                              {
 
  ///---get MDC MC hits and make a map vs wire id
  IDataProviderSvc* eventSvc = NULL;
  Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
  SmartDataPtr<Event::MdcMcHitCol> mdcMcHitCol(eventSvc,"/Event/MC/MdcMcHitCol");
  if(!mdcMcHitCol){
    std::cout<<"MdcUtilitySvc::getMdcMCAssoiciation() does not find MdcMcHitCol!"<<std::endl;
    return;
  }
  std::map<int,Event::MdcMcHit*> mdcMcHitMap;
  Event::MdcMcHitCol::iterator iterMdcMcHit=mdcMcHitCol->begin();
  for(;iterMdcMcHit!=mdcMcHitCol->end();iterMdcMcHit++ )
  {
    int id=m_mdcGeomSvc->Wire(MdcID::layer((*iterMdcMcHit)->identify()),
        MdcID::wire((*iterMdcMcHit)->identify()))->Id();
    mdcMcHitMap.insert(std::pair<int,Event::MdcMcHit*>(id,*iterMdcMcHit));
  }
 
  ///loop over MdcDigi
  MdcDigiVec::const_iterator iterMdcDigi=mdcDigiVecInput.begin();
  for(;iterMdcDigi!=mdcDigiVecInput.end();iterMdcDigi++){
    int digiTrackIndex=(*iterMdcDigi)->getTrackIndex();
    std::cout<<__FILE__<<"   "<<__LINE__<<" "<<digiTrackIndex<<std::endl;
    if(digiTrackIndex>999) digiTrackIndex-=1000;//FIXME
    if(digiTrackIndex!=digiTrackIndex) continue;
    int id=m_mdcGeomSvc->Wire(MdcID::layer((*iterMdcMcHit)->identify()),
        MdcID::wire((*iterMdcMcHit)->identify()))->Id();
    std::map<MdcDigi*, MdcMcHit*> temp;
    std::map<int,Event::MdcMcHit*>::iterator itMdcMcHit=mdcMcHitMap.find(id);
    if(itMdcMcHit!=mdcMcHitMap.end()){
      temp.insert(make_pair(*iterMdcDigi,(*itMdcMcHit).second));
      mdcMCAssociation.insert(make_pair(trackIndex, temp));
      //mdcMCAssociation.insert(std::pair<int,std::map<MdcDigi*,MdcMcHit*> >
      //(trackIndex,make_pair(*iterMdcDigi,mdcMcHitMap[id])));
    }
  }
 
}//end of getMdcMCAssoiciation

getMomPosOfMcParticle()

void MdcUtilitySvc::getMomPosOfMcParticle	(	const Event::McParticle *	mcParticle,
		HepVector3D &	pos,
		HepVector3D &	mom )

Definition at line 910 of file MdcUtilitySvc.cxx.

                                                                                                             {
  if(nullptr==mcParticle) return;
  pos.setX(mcParticle->initialPosition().x());
  pos.setY(mcParticle->initialPosition().y());
  pos.setZ(mcParticle->initialPosition().z());
  mom.setX(mcParticle->initialFourMomentum().px());
  mom.setY(mcParticle->initialFourMomentum().py());
  mom.setZ(mcParticle->initialFourMomentum().pz());
}

Hel()

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

Definition at line 747 of file MdcUtilitySvc.cxx.

                                                                                                                                      {
  double x = piv.x() + dr*cos(phi0) + (Alpha_L/kappa) * (cos(phi0) - cos(phi0+dphi));
  double y = piv.y() + dr*sin(phi0) + (Alpha_L/kappa) * (sin(phi0) - sin(phi0+dphi));
  double z = piv.z() + dz - (Alpha_L/kappa) * dphi * tanl;
  //cout<<"HepPoint3D(x, y, z) = "<<HepPoint3D(x, y, z)<<endl;
  if((x>-1000. && x<1000.) || (y >-1000. && y <1000. ) ||(z>-1000. && z<1000.)){
    return HepPoint3D(x, y, z);
  }else{
    return HepPoint3D(0,0,0);
  }
}

Referenced by cellTrackPassed().

initialize()

StatusCode MdcUtilitySvc::initialize ( )

virtual

Definition at line 32 of file MdcUtilitySvc.cxx.

                                    {
  MsgStream log(messageService(), name());
  log << MSG::INFO << "MdcUtilitySvc::initialize()" << endreq;
 
  StatusCode sc = Service::initialize();
  if( sc.isFailure() ) {
    log << MSG::ERROR << "Service::initialize() failure" << endreq;
    return sc;
  }
 
  //Initalze magnetic flied
 
  sc = service("MagneticFieldSvc", m_pIMF);
  if(! m_pIMF){
    log << MSG::FATAL <<" ERROR Unable to open Magnetic field service "<< endreq;
    return StatusCode::FAILURE;
  }
 
  // Initialize MdcGeomSvc
  IMdcGeomSvc* imdcGeomSvc;
  sc = service("MdcGeomSvc", imdcGeomSvc);
  m_mdcGeomSvc =  dynamic_cast<MdcGeomSvc*> (imdcGeomSvc);
  if(!imdcGeomSvc){
    log << MSG::FATAL <<" Could not load MdcGeomSvc! "<< endreq;
    return StatusCode::FAILURE;
  }
 
  IRawDataProviderSvc* irawDataProviderSvc;
  sc = service ("RawDataProviderSvc", irawDataProviderSvc);
  m_RawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
  if ( sc.isFailure() ){
    log << MSG::FATAL << name()<<" Could not load RawDataProviderSvc!" << endreq;
    return StatusCode::FAILURE;
  }
 
  return StatusCode::SUCCESS;
}

momentum()

Hep3Vector MdcUtilitySvc::momentum

(

const RecMdcTrack *

trk

)

const

Definition at line 825 of file MdcUtilitySvc.cxx.

                                                              {
  // double dr   = trk->helix(0);
  double fi0  = trk->helix(1);
  double cpa  = trk->helix(2);
  double tanl = trk->helix(4);
 
  double pxy = 0.;
  if(cpa != 0) pxy = 1/fabs(cpa);
 
  double px = pxy * (-sin(fi0));
  double py = pxy * cos(fi0);
  double pz = pxy * tanl;
 
  Hep3Vector p;
  p.set(px, py, pz); // MeV/c
  return p;
}

nLayerTrackPassed() [1/2]

int MdcUtilitySvc::nLayerTrackPassed

(

const double

helix[5]

)

const

Definition at line 251 of file MdcUtilitySvc.cxx.

                                                                  {
 
  HepVector helixBesParam(5,0);
  for(int i=0; i<5; ++i) helixBesParam[i] = helixBes[i];
 
  return nLayerTrackPassed(helixBesParam);
}

nLayerTrackPassed() [2/2]

int MdcUtilitySvc::nLayerTrackPassed

(

const HepVector

helix

)

const

Definition at line 260 of file MdcUtilitySvc.cxx.

                                                                  {
  HepVector helix= besPar2PatPar(helixBes);
 
  int nLayer = 0;
 
  for(unsigned iLayer=0; iLayer<43; iLayer++){
    //flightLength is the path length of track in the x-y plane
    //guess flightLength by the radius in middle of the wire.
    double rMidLayer = m_mdcGeomSvc->Layer(iLayer)->Radius();
    double flightLength = rMidLayer;
 
    HepPoint3D pivot(0.,0.,0.);
    double dz = helix[3];
    double c = CLHEP::c_light * 100.; //unit from m/s to cm/s
    double alpha = 1/(c * Bz());//~333.567
    double kappa = helix[2];
    double rc = (-1.)*alpha/kappa;
    //std::cout<<"MdcUtilitySvc alpha   "<<alpha<<std::endl;
    double tanl = helix[4];
    double phi0 = helix[1];
    double phi = flightLength/rc + phi0;//turning angle
    double z = pivot.z() + dz - (alpha/kappa) * tanl * phi;
 
    double layerHalfLength = m_mdcGeomSvc->Layer(iLayer)->Length()/2.;
 
    //std::cout<<"MdcUtilitySvc length  "<<layerHalfLength<<std::endl;
 
    if (fabs(z) < fabs(layerHalfLength)) ++nLayer;
  }
 
  return nLayer;
}

Referenced by nLayerTrackPassed().

p_cms()

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

Definition at line 793 of file MdcUtilitySvc.cxx.

                                                                        {
  HepLorentzVector p4;
  p4.setPx(- sin(helix[1]) / fabs(helix[2]));
  p4.setPy(cos(helix[1]) / fabs(helix[2]));
  p4.setPz(helix[4] / fabs(helix[2]));
  double p3 = p4.mag();
  mass = 0.000511;
  p4.setE(sqrt(p3 * p3 + mass * mass));
 
  double ecm;
  if (runNo > 9815) {
    ecm = 3.097;
  }else{
    ecm = 3.68632;
  }
  double zboost = 0.0075;
  //HepLorentzVector psip(0.011 * 3.68632, 0, 0.0075, 3.68632);
  HepLorentzVector psip(0.011 * ecm, 0, zboost, ecm);
  //cout << setw(15) << ecm << setw(15) << zboost << endl;
  Hep3Vector boostv = psip.boostVector();
 
  //std::cout<<__FILE__<<" boostv "<<boostv<<  std::endl;
  p4.boost(- boostv);
 
  //std::cout<<__FILE__<<" p4 "<<p4<<  std::endl;
  double p_cms = p4.rho();
  //phicms = p4.phi();
  //p4.boost(boostv);
 
  return p_cms;
}

Referenced by p_cms().

patErr2BesErr()

HepSymMatrix MdcUtilitySvc::patErr2BesErr

(

const HepSymMatrix &

err

)

const

Definition at line 311 of file MdcUtilitySvc.cxx.

                                                                      {
  //error matrix
  //std::cout<<" err1  "<<err<<" "<<err.num_row()<<std::endl;
  //V(Y)=S * V(X) * ST , mS = S , mVy = V(Y) , err() = V(X)
  //int n = err.num_row();
  HepSymMatrix mS(err.num_row(),0);
  mS[0][0]=-1.;//dr0=-d0
  mS[1][1]=1.;
 
  mS[2][2]=Bz()/-333.567;//pxy -> cpa
  mS[3][3]=1.;
  mS[4][4]=1.;
  HepSymMatrix mVy= err.similarity(mS);
  //std::cout<<" err2  "<<n<<" "<<mVy<<std::endl;
  return mVy;
}

patPar2BesPar()

HepVector MdcUtilitySvc::patPar2BesPar

(

const HepVector &

helixPar

)

const

Definition at line 294 of file MdcUtilitySvc.cxx.

                                                                     {
  HepVector helix(5,0);
  double d0 = -helixPar[0];    //cm
  double phi0 = helixPar[1]+ CLHEP::halfpi;
  double omega = Bz()*helixPar[2]/-333.567;
  double z0 = helixPar[3];    //cm
  double tanl = helixPar[4];
  helix[0] = d0;
  helix[1] = phi0;
  helix[2] = omega;
  helix[3] = z0;
  helix[4] = tanl;
  //std::cout<<"helix   "<<helix<<std::endl;
  return helix;
}

Referenced by cellTrackPassedPatPar().

pointOnHelix()

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

Definition at line 330 of file MdcUtilitySvc.cxx.

                                                                                                 {
  HepVector helixPat= besPar2PatPar(helixBes);
  return pointOnHelixPatPar(helixPat, layer, innerOrOuter);
}

pointOnHelixPatPar()

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

Definition at line 337 of file MdcUtilitySvc.cxx.

                                                                                                       {
 
  double rInner,rOuter;
  //retrieve Mdc geometry information
  double rCize1 =0.1* m_mdcGeomSvc->Layer(layer)->RCSiz1();  //mm -> cm
  double rCize2 =0.1* m_mdcGeomSvc->Layer(layer)->RCSiz2();  //mm -> cm
  double rLay   =0.1* m_mdcGeomSvc->Layer(layer)->Radius();  //mm -> cm
 
  //(conversion of the units mm(mc)->cm(rec))
  rInner = rLay - rCize1 ; //radius of inner field wire
  rOuter = rLay + rCize2 ; //radius of outer field wire
 
  //int sign = -1; //assumed the first half circle
  HepPoint3D piv(0.,0.,0.);
  double r;
  if (innerOrOuter) r = rInner;
  else r = rOuter;
 
  double d0 = helixPat[0];
  double phi0 = helixPat[1];
  double omega = helixPat[2];
  double z0 = helixPat[3];
  double tanl = helixPat[4];
 
  double rc;
  if (abs(omega) <Constants::epsilon) rc = 9999.;
  else rc = 1./omega;
  double xc     = piv.x() + ( d0 + rc) * cos(phi0);
  double yc     = piv.y() + ( d0 + rc) * sin(phi0);
  HepPoint3D helixCenter(xc,yc,0.0);
  rc = sqrt(xc*xc + yc*yc);//?
  double a,b,c;
  a = r;
  b = d0 + rc;
  c = rc;
  double dphi = acos((a*a-b*b-c*c)/(-2.*b*c));
  double fltlen = dphi * rc;
  double phi = phi0 * omega * fltlen;
  double x = piv.x()+d0*sin(phi) - (rc+d0)*sin(phi0);
  double y = piv.y()+d0*cos(phi) + (rc+d0)*cos(phi0);
  double z = piv.z()+ z0 + fltlen*tanl;
  //std::cout<<" xc yc "<<xc<<" "<<yc
  if(m_debug) {
    std::cout<<" abc "<<a<<" "<<b<<" "<<c<<" omega "<<omega<<" r "<<r<<" rc "<<rc<<" dphi "<<dphi<<" piv  "<<piv.z()
      <<" z0  "<<z0<<" fltlen  "<<fltlen<<" tanl "<<tanl<<std::endl;
    std::cout<<" pointOnHelixPatPar in Hel  "<<helixPat<<std::endl;
    cout<<"HepPoint3D(x, y, z) = "<<HepPoint3D(x, y, z)<<endl;
  }
  return HepPoint3D(x, y, z);
}

Referenced by docaPatPar(), and pointOnHelix().

probab()

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

Definition at line 843 of file MdcUtilitySvc.cxx.

                                                                      {
 
  //constants
  double srtopi=2.0/sqrt(2.0*M_PI);
  double upl=100.0;
 
  double prob=0.0;
  if(ndof<=0) {return prob;}
  if(chisq<0.0) {return prob;}
  if(ndof<=60) {
    //full treatment
    if(chisq>upl) {return prob;}
    double sum=exp(-0.5*chisq);
    double term=sum;
 
    int m=ndof/2;
    if(2*m==ndof){
      if(m==1){return sum;}
      for(int i=2; i<=m;i++){
        term=0.5*term*chisq/(i-1);
        sum+=term;
      }//(int i=2; i<=m)
      return sum;
      //even
 
    }else{
      //odd
      double srty=sqrt(chisq);
      double temp=srty/M_SQRT2;
      prob=erfc(temp);
      if(ndof==1) {return prob;}
      if(ndof==3) {return (srtopi*srty*sum+prob);}
      m=m-1;
      for(int i=1; i<=m; i++){
        term=term*chisq/(2*i+1);
        sum+=term;
      }//(int i=1; i<=m; i++)
      return (srtopi*srty*sum+prob);
 
    }//(2*m==ndof)
 
  }else{
    //asymtotic Gaussian approx
    double srty=sqrt(chisq)-sqrt(ndof-0.5);
    if(srty<12.0) {prob=0.5*erfc(srty);};
 
  }//ndof<30
 
  return prob;
}//endof probab

---

The documentation for this class was generated from the following files:

• MdcUtilitySvc.h
• MdcUtilitySvc.cxx