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 25 of file MdcUtilitySvc.cxx.

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

Referenced by MdcUtilitySvc().

~MdcUtilitySvc()

MdcUtilitySvc::~MdcUtilitySvc

(

)

Definition at line 30 of file MdcUtilitySvc.cxx.

{}

Member Function Documentation

besErr2PatErr()

HepSymMatrix MdcUtilitySvc::besErr2PatErr

(

const HepSymMatrix &

err

)

const

Definition at line 838 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 822 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(), nLayerTrackPassed(), and pointOnHelix().

cellTrackPassed()

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

Definition at line 651 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 590 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 600 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 643 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 401 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(), and docaPatPar().

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 424 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 412 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 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
  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(), and docaPatPar().

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 475 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 456 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 74 of file MdcUtilitySvc.cxx.

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

getChargeOfMcParticle()

float MdcUtilitySvc::getChargeOfMcParticle

(

const Event::McParticle *

mcParticle

)

Definition at line 958 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 982 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 997 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 = nullptr;
      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 81 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 1019 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 971 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 807 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( msgSvc(), name() );
  log << MSG::INFO << "MdcUtilitySvc::initialize()" << endmsg;
 
  StatusCode sc = Service::initialize();
  if ( sc.isFailure() )
  {
    log << MSG::ERROR << "Service::initialize() failure" << endmsg;
    return sc;
  }
 
  // Initalze magnetic flied
 
  sc = service( "MagneticFieldSvc", m_pIMF );
  if ( !m_pIMF )
  {
    log << MSG::FATAL << " ERROR Unable to open Magnetic field service " << endmsg;
    return StatusCode::FAILURE;
  }
 
  // Initialize MdcGeomSvc
  IMdcGeomSvc* imdcGeomSvc;
  sc           = service( "MdcGeomSvc", imdcGeomSvc );
  m_mdcGeomSvc = imdcGeomSvc;
  if ( !imdcGeomSvc )
  {
    log << MSG::FATAL << " Could not load MdcGeomSvc! " << endmsg;
    return StatusCode::FAILURE;
  }
 
  IRawDataProviderSvc* irawDataProviderSvc;
  sc                   = service( "RawDataProviderSvc", irawDataProviderSvc );
  m_RawDataProviderSvc = irawDataProviderSvc;
  if ( sc.isFailure() )
  {
    log << MSG::FATAL << name() << " Could not load RawDataProviderSvc!" << endmsg;
    return StatusCode::FAILURE;
  }
 
  return StatusCode::SUCCESS;
}

momentum()

Hep3Vector MdcUtilitySvc::momentum

(

const RecMdcTrack *

trk

)

const

Definition at line 883 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 257 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 266 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 854 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 318 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 301 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 337 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 345 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 901 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