CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | List of all members
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 ()
 
StatusCode queryInterface (const InterfaceID &riid, void **ppvUnknown)
 
int nLayerTrackPassed (const HepVector helix) const
 
int nLayerTrackPassed (const double helix[5]) const
 
HepVector patPar2BesPar (const HepVector &helixPar) const
 
HepSymMatrix patErr2BesErr (const HepSymMatrix &err) const
 
HepVector besPar2PatPar (const HepVector &helixPar) const
 
HepSymMatrix besErr2PatErr (const HepSymMatrix &err) const
 
double doca (int layer, int cell, const HepVector helix, const HepSymMatrix errMat, bool passCellRequired=true, bool doSag=true) const
 
double doca (int layer, int cell, HepPoint3D eastP, HepPoint3D westP, const HepVector helixBes, const HepSymMatrix errMatBes, bool passCellRequired=true, bool doSag=true) const
 
double doca (int layer, int cell, const MdcSWire *sWire, const HepVector helixPat, const HepSymMatrix errMatPat, bool passCellRequired=true) const
 
double docaPatPar (int layer, int cell, const HepVector helixPat, const HepSymMatrix errMatPat, bool passCellRequired=true, bool doSag=true) const
 
double docaPatPar (int layer, int cell, HepPoint3D eastP, HepPoint3D westP, const HepVector helixBes, const HepSymMatrix errMatBes, bool passCellRequired=true, bool doSag=true) const
 
double docaPatPar (int layer, int cell, const MdcSWire *sWire, const HepVector helixPat, const HepSymMatrix errMatPat, bool passCellRequired=true) const
 
HepPoint3D pointOnHelix (const HepVector helixPar, int lay, int innerOrOuter) const
 
HepPoint3D pointOnHelixPatPar (const HepVector helixPat, int lay, int innerOrOuter) const
 
bool cellTrackPassedByPhi (const HepVector helix, int layer, int &cellId_in, int &cellId_out) const
 
bool cellTrackPassedByPhiPatPar (const HepVector helix, int layer, int &cellId_in, int &cellId_out) const
 
bool cellTrackPassed (const HepVector helix, int layer, int &cellId_in, int &cellId_out) const
 
bool cellTrackPassedPatPar (const HepVector helix, int layer, int &cellId_in, int &cellId_out) const
 
HepPoint3D Hel (HepPoint3D piv, double dr, double phi0, double Alpha_L, double kappa, double dz, double dphi, double tanl) const
 
double p_cms (HepVector helix, int runNo, double mass) const
 
Hep3Vector momentum (const RecMdcTrack *trk) const
 
double probab (const int &ndof, const double &chisq) const
 
virtual int nLayerTrackPassed (const HepVector helix) const =0
 
virtual int nLayerTrackPassed (const double helix[5]) const =0
 
virtual HepVector patPar2BesPar (const HepVector &helixPar) const =0
 
virtual HepSymMatrix patErr2BesErr (const HepSymMatrix &err) const =0
 
virtual HepVector besPar2PatPar (const HepVector &helixPar) const =0
 
virtual HepSymMatrix besErr2PatErr (const HepSymMatrix &err) const =0
 
virtual double doca (int layer, int cell, const HepVector helix, const HepSymMatrix errMat, bool passCellRequired=true, bool doSag=true) const =0
 
virtual double doca (int layer, int cell, HepPoint3D eastP, HepPoint3D westP, const HepVector helixBes, const HepSymMatrix errMatBes, bool passCellRequired=true, bool doSag=true) const =0
 
virtual double doca (int layer, int cell, const MdcSWire *sWire, const HepVector helixPat, const HepSymMatrix errMatPat, bool passCellRequired=true) const =0
 
virtual double docaPatPar (int layer, int cell, const HepVector helixPat, const HepSymMatrix errMatPat, bool passCellRequired=true, bool doSag=true) const =0
 
virtual double docaPatPar (int layer, int cell, HepPoint3D eastP, HepPoint3D westP, const HepVector helixBes, const HepSymMatrix errMatBes, bool passCellRequired=true, bool doSag=true) const =0
 
virtual double docaPatPar (int layer, int cell, const MdcSWire *sWire, const HepVector helixPat, const HepSymMatrix errMatPat, bool passCellRequired=true) const =0
 
virtual HepPoint3D pointOnHelix (const HepVector helixPar, int lay, int innerOrOuter) const =0
 
virtual HepPoint3D pointOnHelixPatPar (const HepVector helixPat, int lay, int innerOrOuter) const =0
 
virtual bool cellTrackPassedByPhi (const HepVector helix, int layer, int &cellId_in, int &cellId_out) const =0
 
virtual bool cellTrackPassedByPhiPatPar (const HepVector helix, int layer, int &cellId_in, int &cellId_out) const =0
 
virtual bool cellTrackPassed (const HepVector helix, int layer, int &cellId_in, int &cellId_out) const =0
 
virtual bool cellTrackPassedPatPar (const HepVector helix, int layer, int &cellId_in, int &cellId_out) const =0
 
virtual HepPoint3D Hel (HepPoint3D piv, double dr, double phi0, double Alpha_L, double kappa, double dz, double dphi, double tanl) const =0
 
virtual double p_cms (HepVector helix, int runNo, double mass) const =0
 
virtual Hep3Vector momentum (const RecMdcTrack *trk) const =0
 
virtual double probab (const int &ndof, const double &chisq) const =0
 

Additional Inherited Members

- Static Public Member Functions inherited from IMdcUtilitySvc
static const InterfaceID & interfaceID ()
 

Detailed Description

Definition at line 19 of file MdcUtilitySvc.h.

Constructor & Destructor Documentation

◆ MdcUtilitySvc()

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

Definition at line 22 of file MdcUtilitySvc.cxx.

22 :
23 Service (name, svcloc) {
24 declareProperty("debug", m_debug = 0);
25 }

◆ ~MdcUtilitySvc()

MdcUtilitySvc::~MdcUtilitySvc ( )

Definition at line 27 of file MdcUtilitySvc.cxx.

27 {
28}

Member Function Documentation

◆ besErr2PatErr()

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

Implements IMdcUtilitySvc.

Definition at line 603 of file MdcUtilitySvc.cxx.

603 {
604 //error matrix
605 //std::cout<<" err1 "<<err<<" "<<err.num_row()<<std::endl;
606 //V(Y)=S * V(X) * ST , mS = S , mVy = V(Y) , err() = V(X)
607 //int n = err.num_row();
608 HepSymMatrix mS(err.num_row(),0);
609 mS[0][0]=-1.;//dr0=-d0
610 mS[1][1]=1.;
611 mS[2][2]=Bz()/-333.567;//pxy -> cpa
612 mS[3][3]=1.;
613 mS[4][4]=1.;
614 HepSymMatrix mVy= err.similarity(mS);
615 //std::cout<<" err2 "<<n<<" "<<mVy<<std::endl;
616 return mVy;
617}

Referenced by doca().

◆ besPar2PatPar()

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

Implements IMdcUtilitySvc.

Definition at line 587 of file MdcUtilitySvc.cxx.

587 {
588 HepVector helix(5,0);
589 double d0 = -helixPar[0]; //cm
590 double phi0 = helixPar[1]+ CLHEP::halfpi;
591 double omega = Bz()*helixPar[2]/-333.567;
592 double z0 = helixPar[3]; //cm
593 double tanl = helixPar[4];
594 helix[0] = d0;
595 helix[1] = phi0;
596 helix[2] = omega;
597 helix[3] = z0;
598 helix[4] = tanl;
599 return helix;
600}

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

◆ cellTrackPassed()

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

Implements IMdcUtilitySvc.

Definition at line 425 of file MdcUtilitySvc.cxx.

425 {
426 int charge,type,nCell;
427 double dr0,phi0,kappa,dz0,tanl;
428 double ALPHA_loc,rho,phi,cosphi0,sinphi0,x_hc,y_hc,z_hc;
429 double dphi0,IO_phi,C_alpha,xx,yy;
430 double inlow,inup,outlow,outup,phi_in,phi_out,phi_bin;
431 double rCize1,rCize2,rLay,length,phioffset,slant,shift;
432 double m_crio[2], phi_io[2], stphi[2],phioff[2],dphi[2];
433
434 dr0 = helixBes[0];
435 phi0 = helixBes[1];
436 kappa = helixBes[2];
437 dz0 = helixBes[3];
438 tanl = helixBes[4];
439
440 ALPHA_loc = 1000/(2.99792458*Bz()); //magnetic field const
441 charge = ( kappa >=0 ) ? 1 : -1 ;
442 rho = ALPHA_loc/kappa ;
443 double pi = Constants::pi;
444 phi = fmod(phi0 + 4*pi , 2*pi);
445 cosphi0 = cos(phi);
446 sinphi0 = (1.0 - cosphi0 ) * (1.0 + cosphi0 );
447 sinphi0 = sqrt(( sinphi0 > 0.) ? sinphi0 : 0.);
448 if( phi > pi ) sinphi0 = -sinphi0 ;
449
450 x_hc = 0. + ( dr0 + rho ) * cosphi0;
451 y_hc = 0. + ( dr0 + rho ) * sinphi0;
452 z_hc = 0. + dz0;
453
454
455 HepPoint3D piv(0.,0.,0.);
456 HepPoint3D hcenter(x_hc,y_hc,0.0);
457
458 double m_c_perp(hcenter.perp());
459 Hep3Vector m_c_unit((HepPoint3D)hcenter.unit());
460 HepPoint3D IO = (-1) * charge * m_c_unit;
461
462 dphi0 = fmod(IO.phi()+4*pi, 2*pi) - phi;
463 IO_phi = fmod(IO.phi()+4*pi, 2*pi);
464
465 if(dphi0 > pi) dphi0 -= 2*pi;
466 else if(dphi0 < -pi) dphi0 += 2*pi;
467
468 phi_io[0] = -(1+charge)*pi/2 - charge*dphi0;
469 phi_io[1] = phi_io[0]+1.5*pi;
470
471
472 bool outFlag = false;
473 //retrieve Mdc geometry information
474 rCize1 = 0.1 * m_mdcGeomSvc->Layer(lay)->RCSiz1(); //mm -> cm
475 rCize2 = 0.1 * m_mdcGeomSvc->Layer(lay)->RCSiz2(); //mm -> cm
476 rLay = 0.1 * m_mdcGeomSvc->Layer(lay)->Radius(); //mm -> cm
477 length = 0.1 * m_mdcGeomSvc->Layer(lay)->Length(); //mm -> cm
478 //double halfLength=length/2.;
479 //(conversion of the units mm(mc)->cm(rec))
480 nCell = m_mdcGeomSvc->Layer(lay)->NCell();
481 phioffset = m_mdcGeomSvc->Layer(lay)->Offset();
482 slant = m_mdcGeomSvc->Layer(lay)->Slant();
483 shift = m_mdcGeomSvc->Layer(lay)->Shift();
484 type = m_mdcGeomSvc->Layer(lay)->Sup()->Type();
485
486 m_crio[0] = rLay - rCize1 ; //radius of inner field wir ( beam pipe)
487 m_crio[1] = rLay + rCize2 ; //radius of outer field wir ( beam pipe)
488
489 int sign = -1; //assumed the first half circle
490
491 Hep3Vector iocand[2];
492 Hep3Vector cell_IO[2];
493
494 for(int ii =0; ii<2; ii++){
495 // By law of cosines to calculate the alpha(up and down field wir_Ge)
496 double cos_alpha = (m_c_perp*m_c_perp + m_crio[ii]*m_crio[ii] - rho*rho)
497 / ( 2 * m_c_perp * m_crio[ii] );
498 if(fabs(cos_alpha)>1&&(ii==0)){
499 cos_alpha = -1;
500 outFlag=true;
501 }
502 if(fabs(cos_alpha)>1&&(ii==1)){
503 cos_alpha = (m_c_perp*m_c_perp + m_crio[0]*m_crio[0] - rho*rho)
504 / ( 2 * m_c_perp * m_crio[0] );
505 C_alpha = 2*pi - acos( cos_alpha);
506 }else {
507 C_alpha = acos( cos_alpha );
508 }
509
510 iocand[ii] = m_c_unit;
511 iocand[ii].rotateZ( charge*sign*C_alpha );
512 iocand[ii] *= m_crio[ii];
513
514 xx = iocand[ii].x() - x_hc ;
515 yy = iocand[ii].y() - y_hc ;
516
517 dphi[ii] = atan2(yy,xx) - phi0 - 0.5*pi*(1-charge);
518 dphi[ii] = fmod( dphi[ii] + 8.0*pi,2*pi);
519
520
521 if( dphi[ii] < phi_io[0] ) {
522 dphi[ii] += 2*pi;
523 }else if( dphi[ii] > phi_io[1] ){ //very very nausea
524 dphi[ii] -= 2*pi;
525 }
526
527 cell_IO[ii] = Hel(piv, dr0, phi, ALPHA_loc, kappa,dz0, dphi[ii], tanl);
528
529 //cout<<" cell_IO["<<ii<<"] : "<<cell_IO[ii]<<endl;
530 if( (cell_IO[ii].x()==0 ) && (cell_IO[ii].y()==0) && (cell_IO[ii].z()==0)) continue;
531 }
532 //if(((fabs(cell_IO[0].z())*10-halfLength)>-7.) && ((fabs(cell_IO[1].z())*10-halfLength)>-7.))return true; //Out sensitive area
533
534 cellId_in = cellId_out = -1 ;
535 phi_in = cell_IO[0].phi();
536 phi_in = fmod ( phi_in + 4 * pi, 2 * pi );
537 phi_out = cell_IO[1].phi();
538 phi_out = fmod ( phi_out + 4 * pi, 2 * pi );
539 phi_bin = 2.0 * pi / nCell ;
540
541 //determine the in/out cell id
542 stphi[0] = shift * phi_bin * (0.5 - cell_IO[0].z()/length);
543 stphi[1] = shift * phi_bin * (0.5 - cell_IO[1].z()/length);
544 //stphi[0],stphi[1] to position fo z axis ,the angle of deflxsion in x_Y
545
546 phioff[0] = phioffset + stphi[0];
547 phioff[1] = phioffset + stphi[1];
548
549 for(int kk = 0; kk<nCell ; kk++){
550 //for stereo layer
551 inlow = phioff[0] + phi_bin*kk - phi_bin*0.5;
552 inlow = fmod( inlow + 4.0 * pi , 2.0 * pi);
553 inup = phioff[0] + phi_bin*kk + phi_bin*0.5;
554 inup = fmod( inup + 4.0 * pi , 2.0 * pi);
555 outlow = phioff[1] + phi_bin*kk - phi_bin*0.5;
556 outlow = fmod( outlow + 4.0 * pi ,2.0 * pi);
557 outup = phioff[1] + phi_bin*kk + phi_bin*0.5;
558 outup = fmod( outup + 4.0 * pi , 2.0 * pi);
559
560 if((phi_in>=inlow && phi_in<=inup)) cellId_in = kk;
561 if((phi_out>=outlow&&phi_out<outup)) cellId_out = kk;
562 if(inlow > inup ){
563 if((phi_in>=inlow&&phi_in<2.0*pi)||(phi_in>=0.0&&phi_in<inup)) cellId_in = kk;
564 }
565 if(outlow>outup){
566 if((phi_out>=outlow&&phi_out<=2.0*pi)||(phi_out>=0.0&&phi_out<outup)) cellId_out = kk;
567 }
568 }//end of nCell loop
569
570 return (cellId_in==cellId_out);
571}
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
length
double length
Definition: CgemLineFit.cxx:106
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
Constants::pi
static const double pi
Definition: Constants.h:38
IMdcGeomSvc::Layer
virtual const MdcGeoLayer *const Layer(unsigned id)=0
MdcGeoLayer::Radius
double Radius(void) const
Definition: MdcGeoLayer.h:160
MdcGeoLayer::Slant
double Slant(void) const
Definition: MdcGeoLayer.h:158
MdcGeoLayer::RCSiz2
double RCSiz2(void) const
Definition: MdcGeoLayer.h:163
MdcGeoLayer::Length
double Length(void) const
Definition: MdcGeoLayer.h:161
MdcGeoLayer::Sup
MdcGeoSuper * Sup(void) const
Definition: MdcGeoLayer.h:174
MdcGeoLayer::Shift
double Shift(void) const
Definition: MdcGeoLayer.h:167
MdcGeoLayer::RCSiz1
double RCSiz1(void) const
Definition: MdcGeoLayer.h:162
MdcGeoLayer::NCell
int NCell(void) const
Definition: MdcGeoLayer.h:165
MdcGeoLayer::Offset
double Offset(void) const
Definition: MdcGeoLayer.h:166
MdcGeoSuper::Type
int Type(void) const
Definition: MdcGeoSuper.h:35
MdcUtilitySvc::Hel
HepPoint3D Hel(HepPoint3D piv, double dr, double phi0, double Alpha_L, double kappa, double dz, double dphi, double tanl) const
Definition: MdcUtilitySvc.cxx:573
pi
const float pi
Definition: vector3.h:133

Referenced by cellTrackPassedPatPar().

◆ cellTrackPassedByPhi()

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

Implements IMdcUtilitySvc.

Definition at line 378 of file MdcUtilitySvc.cxx.

378 {
379 HepVector helixPat = besPar2PatPar(helixBes);
380 return cellTrackPassedByPhiPatPar(helixPat, layer, cellId_in, cellId_out);
381}
MdcUtilitySvc::besPar2PatPar
HepVector besPar2PatPar(const HepVector &helixPar) const
Definition: MdcUtilitySvc.cxx:587
MdcUtilitySvc::cellTrackPassedByPhiPatPar
bool cellTrackPassedByPhiPatPar(const HepVector helix, int layer, int &cellId_in, int &cellId_out) const
Definition: MdcUtilitySvc.cxx:387

◆ cellTrackPassedByPhiPatPar()

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

Implements IMdcUtilitySvc.

Definition at line 387 of file MdcUtilitySvc.cxx.

387 {
388 int nCell = m_mdcGeomSvc->Layer(layer)->NCell();
389
390 double dPhiz = (m_mdcGeomSvc->Wire(layer,0)->Forward().phi() - m_mdcGeomSvc->Wire(layer,0)->Backward().phi())*0.5;
391 double rEnd = m_mdcGeomSvc->Wire(layer,0)->Backward().rho()/10.;//mm2cm
392 double rMid = rEnd * cos(dPhiz);
393 //std::cout<<"( cell "<<0<<" westPphi "<<m_mdcGeomSvc->Wire(layer,0)->Forward().phi() <<" eastPphi "
394 //<<m_mdcGeomSvc->Wire(layer,0)->Backward().phi()<<" twist "<<dPhiz<<" rend "<<rEnd<<std::endl;
395 double fltLenIn = rMid;
396 double phiIn = helixPat[1] + helixPat[2]*fltLenIn;//phi0 + omega * L
397
398 double phiEPOffset = m_mdcGeomSvc->Wire(layer,0)->Backward().phi();//east.phi()= BackWard.phi
399 BesAngle tmp(phiIn - phiEPOffset);
400 if(m_debug) std::cout<<"cellTrackPassed nCell "<<nCell<<" layer "<<layer<<" fltLenIn "<<fltLenIn<<" phiEPOffset "<<phiEPOffset<<std::endl;
401 //BesAngle tmp(phiIn - layPtr->phiEPOffset());
402 int wlow = (int)floor(nCell * tmp.rad() / twopi );
403 int wbig = (int)ceil(nCell * tmp.rad() / twopi );
404 if (tmp == 0 ){ wlow = -1; wbig = 1; }
405
406 if ((wlow%nCell)< 0) wlow += nCell;
407 cellId_in = wlow;
408
409 if ((wbig%nCell)< 0) wbig += nCell;
410 cellId_out = wbig;
411
412 bool passedOneCell = (cellId_in == cellId_out);
413
414 return passedOneCell;//pass more than one cell
415}
IMdcGeomSvc::Wire
virtual const MdcGeoWire *const Wire(unsigned id)=0
MdcGeoWire::Forward
HepPoint3D Forward(void) const
Definition: MdcGeoWire.h:129
MdcGeoWire::Backward
HepPoint3D Backward(void) const
Definition: MdcGeoWire.h:128

Referenced by cellTrackPassedByPhi(), and docaPatPar().

◆ cellTrackPassedPatPar()

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

Implements IMdcUtilitySvc.

Definition at line 418 of file MdcUtilitySvc.cxx.

418 {
419 HepVector helixBes = patPar2BesPar(helixPat);
420 return cellTrackPassed(helixBes, layer, cellId_in, cellId_out);
421}
MdcUtilitySvc::patPar2BesPar
HepVector patPar2BesPar(const HepVector &helixPar) const
Definition: MdcUtilitySvc.cxx:120
MdcUtilitySvc::cellTrackPassed
bool cellTrackPassed(const HepVector helix, int layer, int &cellId_in, int &cellId_out) const
Definition: MdcUtilitySvc.cxx:425

◆ doca() [1/3]

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

Implements IMdcUtilitySvc.

Definition at line 216 of file MdcUtilitySvc.cxx.

216 {
217 HepVector helixPat = besPar2PatPar(helixBes);
218 HepSymMatrix errMatPat= besErr2PatErr(errMatBes);
219
220 return docaPatPar(layer, cell, helixPat, errMatPat, passCellRequired, doSag);//call 4.
221}
MdcUtilitySvc::besErr2PatErr
HepSymMatrix besErr2PatErr(const HepSymMatrix &err) const
Definition: MdcUtilitySvc.cxx:603
MdcUtilitySvc::docaPatPar
double docaPatPar(int layer, int cell, const HepVector helixPat, const HepSymMatrix errMatPat, bool passCellRequired=true, bool doSag=true) const
Definition: MdcUtilitySvc.cxx:243

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

◆ doca() [2/3]

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

Implements IMdcUtilitySvc.

Definition at line 234 of file MdcUtilitySvc.cxx.

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

◆ 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
virtual

Implements IMdcUtilitySvc.

Definition at line 225 of file MdcUtilitySvc.cxx.

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

◆ docaPatPar() [1/3]

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

Implements IMdcUtilitySvc.

Definition at line 243 of file MdcUtilitySvc.cxx.

243 {
244
245 const MdcGeoWire *geowir = m_mdcGeomSvc->Wire(layer,cell);
246 int id = geowir->Id();
247 double sag = 0.;
248 if(doSag) sag = geowir->Sag()/10.;//mm2cm
249 HepPoint3D eastP = geowir->Backward()/10.0;//mm2cm
250 HepPoint3D westP = geowir->Forward() /10.0;//mm2cm
251 const MdcSWire* sWire = new MdcSWire(eastP, westP, sag, id , cell);
252
253 double doca = docaPatPar(layer, cell, sWire, helixPat, errMat, passCellRequired);//call 6.
254
255 delete sWire;
256 return doca;
257}
MdcGeoWire::Sag
const double Sag(void) const
Definition: MdcGeoWire.cxx:12
MdcGeoWire::Id
int Id(void) const
Definition: MdcGeoWire.h:127
MdcUtilitySvc::doca
double doca(int layer, int cell, const HepVector helix, const HepSymMatrix errMat, bool passCellRequired=true, bool doSag=true) const
Definition: MdcUtilitySvc.cxx:216

Referenced by doca(), 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
virtual

Implements IMdcUtilitySvc.

Definition at line 277 of file MdcUtilitySvc.cxx.

277 {
278
279 if(m_debug) std::cout<<" helixPat "<<helixPat<<std::endl;
280 if(m_debug) std::cout<<" layer "<<layer<<" cell "<<cell<<std::endl;
281 //-----cell track passed
282 int cellId_in = -1;
283 int cellId_out = -1;
284
285 //cellTrackPassedPatPar(helixPat,layer,cellId_in,cellId_out);//yzhang FIXME 2012-07-11
286 cellTrackPassedByPhiPatPar(helixPat,layer,cellId_in,cellId_out);//yzhang FIXME 2012-07-11
287
288 if(m_debug) {
289 std::cout<<" cellId in "<<cellId_in<<" out "<<cellId_out <<std::endl;
290 cout << "cell = " << cell << ", cellId_in = " << cellId_in << ", cellId_out = " << cellId_out << endl;
291 }
292 if (passCellRequired &&(cell < cellId_in && cell > cellId_out)) return -999.;
293
294 //-----helix trajectory
295 HelixTraj* helixTraj = new HelixTraj(helixPat,errMat);
296 const Trajectory* trajHelix = dynamic_cast<const Trajectory*> (helixTraj);
297
298 //-----pointOnHelix
299 int innerOrOuter = 1;
300 Hep3Vector cell_IO = pointOnHelixPatPar(helixPat,layer,innerOrOuter);
301 double fltTrack = 0.1 * m_mdcGeomSvc -> Layer(layer)->Radius();
302 if(m_debug) {
303 std::cout<<" cell_IO "<<cell_IO<<std::endl;
304 std::cout<<" fltTrack "<<fltTrack<<std::endl;
305 }
306
307 //------wire trajectory
308 const MdcSagTraj* m_wireTraj = sWire->getTraj();
309 const Trajectory* trajWire = dynamic_cast<const Trajectory*> (m_wireTraj);
310 const HepPoint3D* start_In = sWire->getEastPoint();
311 //const HepPoint3D* stop_In = sWire->getWestPoint();
312 if(m_debug){
313 std::cout<<" sag "<<m_wireTraj->sag()<<std::endl;
314 std::cout<< " east -------- "<< start_In->x()<<","<<start_In->y()<<","<<start_In->z()<<std::endl;
315 }
316 //std::cout<< " west -------- "<< stop_In->x()<<","<<stop_In->y()<<","<<stop_In->z() <<std::endl;
317
318 //------calc poca
319 double doca = -999.;
320 TrkPoca* trkPoca;
321 double zWire = cell_IO.z();
322 //HepPoint3D pos_in(zWire,sWire->xWireDC(zWire),sWire->yWireDC(zWire)) ;
323 HepPoint3D pos_in(sWire->xWireDC(zWire),sWire->yWireDC(zWire),zWire) ;
324 if(m_debug) std::cout<<" zWire "<<zWire<<" zEndDC "<<sWire->zEndDC()<<std::endl;
325 //if(m_debug) std::cout<<"pos_in "<<pos_in<<" fltWire "<<fltWire<<std::endl;
326
327 double fltWire = sqrt( (pos_in.x()-start_In->x())*(pos_in.x()-start_In->x()) +
328 (pos_in.y()-start_In->y())*(pos_in.y()-start_In->y()) +
329 (pos_in.z()-start_In->z())*(pos_in.z()-start_In->z()) );
330 trkPoca = new TrkPoca(*trajHelix, fltTrack, *trajWire, fltWire);
331 doca = trkPoca->doca();
332
333 double hitLen = trkPoca->flt2();
334 double startLen = trajWire->lowRange() - 5.;
335 double endLen = trajWire->hiRange() + 5.;
336 if(hitLen < startLen || hitLen > endLen) {
337 if(m_debug) std::cout<<"WARNING MdcUtilitySvc::docaPatPar() isBeyondEndflange! hitLen="<<hitLen <<" startLen="<<startLen<<" endLen "<<endLen<<std::endl;
338 doca = -998.;
339 }
340 //std::cout<<" hitLen "<<hitLen<<" startLen "<<startLen<<" endLen "<<endLen <<" doca "<<doca<<std::endl;
341
342 if(m_debug) std::cout<<" doca "<<doca<<std::endl;
343 delete helixTraj;
344 delete trkPoca;
345
346 return doca;
347
348} //----------end of calculatng the doca ---------------------------------//
MdcSWire::yWireDC
double yWireDC(double z) const
Definition: MdcSWire.h:62
MdcSWire::xWireDC
double xWireDC(double z) const
Definition: MdcSWire.h:61
MdcSWire::zEndDC
double zEndDC(void) const
Definition: MdcSWire.h:49
MdcSWire::getTraj
const MdcSagTraj * getTraj(void) const
Definition: MdcSWire.h:34
MdcSWire::getEastPoint
const HepPoint3D * getEastPoint(void) const
Definition: MdcSWire.h:32
MdcSagTraj::sag
double sag(void) const
Definition: MdcSagTraj.h:80
MdcUtilitySvc::pointOnHelixPatPar
HepPoint3D pointOnHelixPatPar(const HepVector helixPat, int lay, int innerOrOuter) const
Definition: MdcUtilitySvc.cxx:163
Trajectory::lowRange
double lowRange() const
Definition: Trajectory.h:91
Trajectory::hiRange
double hiRange() const
Definition: Trajectory.h:92
TrkPocaBase::flt2
double flt2() const
Definition: TrkPocaBase.h:68
TrkPoca::doca
double doca() const
Definition: TrkPoca.h:56

◆ 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
virtual

Implements IMdcUtilitySvc.

Definition at line 261 of file MdcUtilitySvc.cxx.

261 {
262
263 const MdcGeoWire *geowir = m_mdcGeomSvc->Wire(layer,cell);
264 int id = geowir->Id();
265 double sag = 0.;
266 if(doSag) sag = geowir->Sag()/10.;//mm2cm
267 const MdcSWire* sWire = new MdcSWire(eastP, westP, sag, id , cell);//cm
268
269 double doca = docaPatPar(layer, cell, sWire, helixPat, errMat, passCellRequired);//call 6.
270
271 delete sWire;
272 return doca;
273}

◆ finalize()

StatusCode MdcUtilitySvc::finalize ( )
virtual

Definition at line 57 of file MdcUtilitySvc.cxx.

57 {
58 MsgStream log(messageService(), name());
59 log << MSG::INFO << "MdcUtilitySvc::finalize()" << endreq;
60
61 return StatusCode::SUCCESS;
62}

◆ Hel()

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

Implements IMdcUtilitySvc.

Definition at line 573 of file MdcUtilitySvc.cxx.

573 {
574 double x = piv.x() + dr*cos(phi0) + (Alpha_L/kappa) * (cos(phi0) - cos(phi0+dphi));
575 double y = piv.y() + dr*sin(phi0) + (Alpha_L/kappa) * (sin(phi0) - sin(phi0+dphi));
576 double z = piv.z() + dz - (Alpha_L/kappa) * dphi * tanl;
577 //cout<<"HepPoint3D(x, y, z) = "<<HepPoint3D(x, y, z)<<endl;
578 if((x>-1000. && x<1000.) || (y >-1000. && y <1000. ) ||(z>-1000. && z<1000.)){
579 return HepPoint3D(x, y, z);
580 }else{
581 return HepPoint3D(0,0,0);
582 }
583}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: Gam4pikp.cxx:37

Referenced by cellTrackPassed().

◆ initialize()

StatusCode MdcUtilitySvc::initialize ( )
virtual

Definition at line 30 of file MdcUtilitySvc.cxx.

30 {
31 MsgStream log(messageService(), name());
32 log << MSG::INFO << "MdcUtilitySvc::initialize()" << endreq;
33
34 StatusCode sc = Service::initialize();
35 if( sc.isFailure() ) {
36 log << MSG::ERROR << "Service::initialize() failure" << endreq;
37 return sc;
38 }
39
40 //Initalze magnetic flied
41 sc = service("MagneticFieldSvc", m_pIMF);
42 if(! m_pIMF){
43 log << MSG::FATAL <<" ERROR Unable to open Magnetic field service "<< endreq;
44 return StatusCode::FAILURE;
45 }
46
47 // Initialize MdcGeomSvc
48 sc = service("MdcGeomSvc", m_mdcGeomSvc);
49 if(! m_mdcGeomSvc){
50 log << MSG::FATAL <<" Could not load MdcGeomSvc! "<< endreq;
51 return StatusCode::FAILURE;
52 }
53
54 return StatusCode::SUCCESS;
55}

◆ momentum()

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

Implements IMdcUtilitySvc.

Definition at line 651 of file MdcUtilitySvc.cxx.

651 {
652 // double dr = trk->helix(0);
653 double fi0 = trk->helix(1);
654 double cpa = trk->helix(2);
655 double tanl = trk->helix(4);
656
657 double pxy = 0.;
658 if(cpa != 0) pxy = 1/fabs(cpa);
659
660 double px = pxy * (-sin(fi0));
661 double py = pxy * cos(fi0);
662 double pz = pxy * tanl;
663
664 Hep3Vector p;
665 p.set(px, py, pz); // MeV/c
666 return p;
667}
DstMdcTrack::helix
const HepVector helix() const
......
Definition: DstMdcTrack.cxx:157

◆ nLayerTrackPassed() [1/2]

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

Implements IMdcUtilitySvc.

Definition at line 77 of file MdcUtilitySvc.cxx.

77 {
78
79 HepVector helixBesParam(5,0);
80 for(int i=0; i<5; ++i) helixBesParam[i] = helixBes[i];
81
82 return nLayerTrackPassed(helixBesParam);
83}
MdcUtilitySvc::nLayerTrackPassed
int nLayerTrackPassed(const HepVector helix) const
Definition: MdcUtilitySvc.cxx:86

◆ nLayerTrackPassed() [2/2]

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

Implements IMdcUtilitySvc.

Definition at line 86 of file MdcUtilitySvc.cxx.

86 {
87 HepVector helix= besPar2PatPar(helixBes);
88
89 int nLayer = 0;
90
91 for(unsigned iLayer=0; iLayer<43; iLayer++){
92 //flightLength is the path length of track in the x-y plane
93 //guess flightLength by the radius in middle of the wire.
94 double rMidLayer = m_mdcGeomSvc->Layer(iLayer)->Radius();
95 double flightLength = rMidLayer;
96
97 HepPoint3D pivot(0.,0.,0.);
98 double dz = helix[3];
99 double c = CLHEP::c_light * 100.; //unit from m/s to cm/s
100 double alpha = 1/(c * Bz());//~333.567
101 double kappa = helix[2];
102 double rc = (-1.)*alpha/kappa;
103 //std::cout<<"MdcUtilitySvc alpha "<<alpha<<std::endl;
104 double tanl = helix[4];
105 double phi0 = helix[1];
106 double phi = flightLength/rc + phi0;//turning angle
107 double z = pivot.z() + dz - (alpha/kappa) * tanl * phi;
108
109 double layerHalfLength = m_mdcGeomSvc->Layer(iLayer)->Length()/2.;
110
111 //std::cout<<"MdcUtilitySvc length "<<layerHalfLength<<std::endl;
112
113 if (fabs(z) < fabs(layerHalfLength)) ++nLayer;
114 }
115
116 return nLayer;
117}
alpha
const double alpha
Definition: FastVertexFit.cxx:4

Referenced by nLayerTrackPassed().

◆ p_cms()

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

Implements IMdcUtilitySvc.

Definition at line 619 of file MdcUtilitySvc.cxx.

619 {
620 HepLorentzVector p4;
621 p4.setPx(- sin(helix[1]) / fabs(helix[2]));
622 p4.setPy(cos(helix[1]) / fabs(helix[2]));
623 p4.setPz(helix[4] / fabs(helix[2]));
624 double p3 = p4.mag();
625 mass = 0.000511;
626 p4.setE(sqrt(p3 * p3 + mass * mass));
627
628 double ecm;
629 if (runNo > 9815) {
630 ecm = 3.097;
631 }else{
632 ecm = 3.68632;
633 }
634 double zboost = 0.0075;
635 //HepLorentzVector psip(0.011 * 3.68632, 0, 0.0075, 3.68632);
636 HepLorentzVector psip(0.011 * ecm, 0, zboost, ecm);
637 //cout << setw(15) << ecm << setw(15) << zboost << endl;
638 Hep3Vector boostv = psip.boostVector();
639
640 //std::cout<<__FILE__<<" boostv "<<boostv<< std::endl;
641 p4.boost(- boostv);
642
643 //std::cout<<__FILE__<<" p4 "<<p4<< std::endl;
644 double p_cms = p4.rho();
645 //phicms = p4.phi();
646 //p4.boost(boostv);
647
648 return p_cms;
649}
mass
double mass
Definition: CosmicGenerator.cxx:138
runNo
int runNo
Definition: DedxCalibRunByRun.cxx:14
MdcUtilitySvc::p_cms
double p_cms(HepVector helix, int runNo, double mass) const
Definition: MdcUtilitySvc.cxx:619

Referenced by p_cms().

◆ patErr2BesErr()

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

Implements IMdcUtilitySvc.

Definition at line 137 of file MdcUtilitySvc.cxx.

137 {
138 //error matrix
139 //std::cout<<" err1 "<<err<<" "<<err.num_row()<<std::endl;
140 //V(Y)=S * V(X) * ST , mS = S , mVy = V(Y) , err() = V(X)
141 //int n = err.num_row();
142 HepSymMatrix mS(err.num_row(),0);
143 mS[0][0]=-1.;//dr0=-d0
144 mS[1][1]=1.;
145
146 mS[2][2]=Bz()/-333.567;//pxy -> cpa
147 mS[3][3]=1.;
148 mS[4][4]=1.;
149 HepSymMatrix mVy= err.similarity(mS);
150 //std::cout<<" err2 "<<n<<" "<<mVy<<std::endl;
151 return mVy;
152}

◆ patPar2BesPar()

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

Implements IMdcUtilitySvc.

Definition at line 120 of file MdcUtilitySvc.cxx.

120 {
121 HepVector helix(5,0);
122 double d0 = -helixPar[0]; //cm
123 double phi0 = helixPar[1]+ CLHEP::halfpi;
124 double omega = Bz()*helixPar[2]/-333.567;
125 double z0 = helixPar[3]; //cm
126 double tanl = helixPar[4];
127 helix[0] = d0;
128 helix[1] = phi0;
129 helix[2] = omega;
130 helix[3] = z0;
131 helix[4] = tanl;
132 //std::cout<<"helix "<<helix<<std::endl;
133 return helix;
134}

Referenced by cellTrackPassedPatPar().

◆ pointOnHelix()

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

Implements IMdcUtilitySvc.

Definition at line 156 of file MdcUtilitySvc.cxx.

156 {
157 HepVector helixPat= besPar2PatPar(helixBes);
158 return pointOnHelixPatPar(helixPat, layer, innerOrOuter);
159}

◆ pointOnHelixPatPar()

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

Implements IMdcUtilitySvc.

Definition at line 163 of file MdcUtilitySvc.cxx.

163 {
164
165 double rInner,rOuter;
166 //retrieve Mdc geometry information
167 double rCize1 =0.1* m_mdcGeomSvc->Layer(layer)->RCSiz1(); //mm -> cm
168 double rCize2 =0.1* m_mdcGeomSvc->Layer(layer)->RCSiz2(); //mm -> cm
169 double rLay =0.1* m_mdcGeomSvc->Layer(layer)->Radius(); //mm -> cm
170
171 //(conversion of the units mm(mc)->cm(rec))
172 rInner = rLay - rCize1 ; //radius of inner field wire
173 rOuter = rLay + rCize2 ; //radius of outer field wire
174
175 //int sign = -1; //assumed the first half circle
176 HepPoint3D piv(0.,0.,0.);
177 double r;
178 if (innerOrOuter) r = rInner;
179 else r = rOuter;
180
181 double d0 = helixPat[0];
182 double phi0 = helixPat[1];
183 double omega = helixPat[2];
184 double z0 = helixPat[3];
185 double tanl = helixPat[4];
186
187 double rc;
188 if (abs(omega) <Constants::epsilon) rc = 9999.;
189 else rc = 1./omega;
190 double xc = piv.x() + ( d0 + rc) * cos(phi0);
191 double yc = piv.y() + ( d0 + rc) * sin(phi0);
192 HepPoint3D helixCenter(xc,yc,0.0);
193 rc = sqrt(xc*xc + yc*yc);//?
194 double a,b,c;
195 a = r;
196 b = d0 + rc;
197 c = rc;
198 double dphi = acos((a*a-b*b-c*c)/(-2.*b*c));
199 double fltlen = dphi * rc;
200 double phi = phi0 * omega * fltlen;
201 double x = piv.x()+d0*sin(phi) - (rc+d0)*sin(phi0);
202 double y = piv.y()+d0*cos(phi) + (rc+d0)*cos(phi0);
203 double z = piv.z()+ z0 + fltlen*tanl;
204 //std::cout<<" xc yc "<<xc<<" "<<yc
205 if(m_debug) {
206 std::cout<<" abc "<<a<<" "<<b<<" "<<c<<" omega "<<omega<<" r "<<r<<" rc "<<rc<<" dphi "<<dphi<<" piv "<<piv.z()
207 <<" z0 "<<z0<<" fltlen "<<fltlen<<" tanl "<<tanl<<std::endl;
208 std::cout<<" pointOnHelixPatPar in Hel "<<helixPat<<std::endl;
209 cout<<"HepPoint3D(x, y, z) = "<<HepPoint3D(x, y, z)<<endl;
210 }
211 return HepPoint3D(x, y, z);
212}
abs
double abs(const EvtComplex &c)
Definition: EvtComplex.hh:212
Constants::epsilon
static const double epsilon
Definition: Constants.h:46

Referenced by docaPatPar(), and pointOnHelix().

◆ probab()

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

Implements IMdcUtilitySvc.

Definition at line 669 of file MdcUtilitySvc.cxx.

669 {
670
671 //constants
672 double srtopi=2.0/sqrt(2.0*M_PI);
673 double upl=100.0;
674
675 double prob=0.0;
676 if(ndof<=0) {return prob;}
677 if(chisq<0.0) {return prob;}
678 if(ndof<=60) {
679 //full treatment
680 if(chisq>upl) {return prob;}
681 double sum=exp(-0.5*chisq);
682 double term=sum;
683
684 int m=ndof/2;
685 if(2*m==ndof){
686 if(m==1){return sum;}
687 for(int i=2; i<=m;i++){
688 term=0.5*term*chisq/(i-1);
689 sum+=term;
690 }//(int i=2; i<=m)
691 return sum;
692 //even
693
694 }else{
695 //odd
696 double srty=sqrt(chisq);
697 double temp=srty/M_SQRT2;
698 prob=erfc(temp);
699 if(ndof==1) {return prob;}
700 if(ndof==3) {return (srtopi*srty*sum+prob);}
701 m=m-1;
702 for(int i=1; i<=m; i++){
703 term=term*chisq/(2*i+1);
704 sum+=term;
705 }//(int i=1; i<=m; i++)
706 return (srtopi*srty*sum+prob);
707
708 }//(2*m==ndof)
709
710 }else{
711 //asymtotic Gaussian approx
712 double srty=sqrt(chisq)-sqrt(ndof-0.5);
713 if(srty<12.0) {prob=0.5*erfc(srty);};
714
715 }//ndof<30
716
717 return prob;
718}//endof probab
exp
EvtComplex exp(const EvtComplex &c)
Definition: EvtComplex.hh:252
M_PI
#define M_PI
Definition: TConstant.h:4

◆ queryInterface()

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

Definition at line 64 of file MdcUtilitySvc.cxx.

64 {
65
66 if( IID_IMdcUtilitySvc.versionMatch(riid) ){
67 *ppvInterface = static_cast<IMdcUtilitySvc*> (this);
68 } else{
69 return Service::queryInterface(riid, ppvInterface);
70 }
71
72 addRef();
73 return StatusCode::SUCCESS;
74}
The documentation for this class was generated from the following files: