32 const std::string &
name,
36 declareInterface<IDerivCalcTool>(
this);
44 return StatusCode::FAILURE;
50 return StatusCode::FAILURE;
54 return StatusCode::SUCCESS;
61 return StatusCode::SUCCESS;
73 std::vector<bool> hitModules(
nModules,
false);
77 std::vector<AlignModule *> alignModules;
79 for (; atsosItr != alignTrack->
lastAtsos(); ++atsosItr) {
86 if (!(*atsosItr)->isValid() || !
module)
90 if(!hitModules[
module->identifyHash()]) {
91 hitModules[
module->identifyHash()] =
true;
92 alignModules.push_back(
module);
97 std::vector<AlignModuleDerivatives> * derivatives =
new std::vector<AlignModuleDerivatives>;
99 for ( ; moduleIt!=alignModules.end(); ++moduleIt) {
100 std::vector<Amg::VectorX> deriv_vec =
getDerivatives(alignTrack,*moduleIt);
101 derivatives->push_back(make_pair(*moduleIt,deriv_vec));
116 static std::once_flag
flag;
117 std::call_once(
flag, [&]() {
119 *
m_logStream<<
"*************************************************************"<<std::endl;
120 *
m_logStream<<
"*************************************************************"<<std::endl;
124 *
m_logStream<<
"*** Running full LOCAL Chi2 method *****"<<std::endl;
126 *
m_logStream<<
"*** Running full GLOBAL Chi2 method *****"<<std::endl;
128 *
m_logStream<<
"*** using analytical derivatives *****"<<std::endl;
131 std::string resType =
"HitOnly ";
133 resType =
"Unbiased";
134 *
m_logStream<<
"*** with residual type: "<<resType<<
" *****"<<std::endl;
138 *
m_logStream<<
"*************************************************************"<<std::endl;
145 ATH_MSG_WARNING(
"Inverse local error matrix is invalid, skipping the track");
150 for (
int i=0;
i<Vinv->rows();
i++)
msg()<<(*Vinv)(
i,
i)<<
" ";
158 ATH_MSG_DEBUG(
"setting Residual covariance matrix for local method");
169 ATH_MSG_DEBUG(
"Ignoring measured errors for Pixel clusters, using intrinsic errors");
171 ATH_MSG_DEBUG(
"Ignoring measured errors for SCT clusters, using intrinsic errors");
173 ATH_MSG_DEBUG(
"Ignoring measured errors for TRT clusters, using intrinsic errors");
178 for ( ; itAtsos != itAtsos_end; ++itAtsos) {
196 ATH_MSG_DEBUG(
"setting local weight matrix W ( "<<
W->rows()<<
" x "<<
W->cols()<<
" ) (diagonal only):");
198 for (
int i=0;
i<
W->rows();
i++)
msg()<<(*W)(
i,
i)<<
" ";
215 ATH_MSG_DEBUG(
"V ( "<<V->rows()<<
" x "<<V->cols()<<
" ) (diagonal only):");
217 for (
int i=0;
i<V->rows();
i++)
msg()<<(*V)(
i,
i)<<
" ";
240 *
W = ((*Vinv) * R * (*Vinv));
243 ATH_MSG_DEBUG(
"Weight matrix is invalid, skipping the track");
253 ATH_MSG_DEBUG(
"setting weight for 1st derivatives (diagonal only): ");
255 for (
int i=0;
i<W1st->rows();
i++)
msg()<<(*W1st)(
i,
i)<<
" ";
276 if( H0==
nullptr || C==
nullptr) {
277 ATH_MSG_ERROR(
"no derivative matrix or cov matrix stored on AlignTrack!"
278 <<
"This should have been done in AlignTrackPreProcessor!"
289 int Csize(C->rows());
304 for(
int ii=0; ii<Csize; ++ii ) {
305 for(
int jj=ii; jj<Csize; ++jj ) {
306 CA->elemr(ii,jj) = (*C)(ii,jj);
313 ATH_MSG_ERROR(
"First inversion of matrix CA failed with LAPACK status flag " << ierr);
318 if( !(alignTrack->
refitD0()) ) {
319 for(
int ii=0; ii<(Csize); ++ii)
CA->elemr(0,ii)=0.0;
322 if( !(alignTrack->
refitZ0()) ) {
323 for(
int ii=0; ii<(Csize); ++ii)
CA->elemr(1,ii)=0.0;
327 for(
int ii=0; ii<(Csize); ++ii)
CA->elemr(2,ii)=0.0;
331 for(
int ii=0; ii<(Csize); ++ii)
CA->elemr(3,ii)=0.0;
335 for(
int ii=0; ii<(Csize); ++ii)
CA->elemr(4,ii)=0.0;
343 ATH_MSG_ERROR(
"Second inversion of matrix CA failed with LAPACK status flag " << ierr);
348 for(
int ii=0; ii<Csize; ++ii ) {
349 for(
int jj=ii; jj<Csize; ++jj ) {
350 CC(ii,jj) =
CA->elemc(ii,jj);
355 if( !(alignTrack->
refitD0()) )
for(
int ii=0; ii<(Csize); ++ii){
CC(0,ii)=0.0;
CC(ii,0)=0.0; };
356 if( !(alignTrack->
refitZ0()) )
for(
int ii=0; ii<(Csize); ++ii){
CC(1,ii)=0.0;
CC(ii,1)=0.0; };
357 if( !(alignTrack->
refitPhi()) )
for(
int ii=0; ii<(Csize); ++ii){
CC(2,ii)=0.0;
CC(ii,2)=0.0; };
358 if( !(alignTrack->
refitTheta()) )
for(
int ii=0; ii<(Csize); ++ii){
CC(3,ii)=0.0;
CC(ii,3)=0.0; };
359 if( !(alignTrack->
refitQovP()) )
for(
int ii=0; ii<(Csize); ++ii){
CC(4,ii)=0.0;
CC(ii,4)=0.0; };
370 int nMeas = H0->rows();
376 HCH =
CC.similarity( *H0 );
378 ATH_MSG_DEBUG(
"HCH ( "<<HCH.rows()<<
" x "<<HCH.cols()<<
" )");
385 std::vector<int> matrixIndices(nMeas);
410 matrixIndices[imeas]=-1;
417 for (
const AlignTSOS* atsos : *alignTSOSCollection) {
419 if (!atsos->isValid())
423 ATH_MSG_ERROR(
"can't use scatterers on AlignTrack yet for analytical derivatives!");
429 if (atsos_rio->
identify()==tsosId) {
430 matrixIndices[imeas]=iameas;
434 if (atsos->nResDim()>1)
438 matrixIndices[imeas]=iameas;
444 else if(atsos->nResDim()>1)
454 matrixIndices[imeas]=-1;
455 ATH_MSG_DEBUG(
"matrixIndices["<<imeas<<
"]="<<matrixIndices[imeas]);
465 for (
int k=0;
k<nMeas;
k++) {
467 int iameas=matrixIndices[
k];
471 for (
int l=0;
l<nMeas;
l++) {
472 int jameas=matrixIndices[
l];
476 Q(iameas,jameas) = HCH(
k,
l);
481 ATH_MSG_DEBUG(
"before check Q ( "<<Q.rows()<<
" x "<<Q.cols()<<
" )");
485 ATH_MSG_DEBUG(
"Matrix Q = HCH is invalid, skipping the track");
498 for (; itAtsos != alignTrack->
lastAtsos(); ++itAtsos) {
500 std::vector<Residual>::const_iterator itRes=(**itAtsos).firstResidual();
501 for (; itRes!=(**itAtsos).lastResidual(); ++itRes,
index++) {
519 const double epsilon=1
e-10;
520 for(
int irow=0; irow<R.rows(); ++irow) {
522 Risvalid = Risvalid && R(irow,irow)>0;
524 if( !(R(irow,irow)>0) )
525 msg(
MSG::DEBUG) <<
"matrix invalid: (" << irow <<
"," << irow<<
") = " << R(irow,irow) <<
endmsg;
532 double Rcorr = R(irow,
icol)/sqrt(R(irow,irow)*R(
icol,
icol));
533 if( Rcorr+epsilon<-1 || Rcorr-epsilon>1 )
537 ATH_MSG_DEBUG(
"matrix corr invalid for (" << irow <<
"," <<
icol <<
") Rcorr = " << Rcorr);
558 globalFrameToAlignFrame(0,1)<<
" "<<
559 globalFrameToAlignFrame(0,2));
561 globalFrameToAlignFrame(1,1)<<
" "<<
562 globalFrameToAlignFrame(1,2));
564 globalFrameToAlignFrame(2,1)<<
" "<<
565 globalFrameToAlignFrame(2,2));
570 globalToAlignFrameRotation(0,1)<<
" "<<
571 globalToAlignFrameRotation(0,2));
573 globalToAlignFrameRotation(1,1)<<
" "<<
574 globalToAlignFrameRotation(1,2));
576 globalToAlignFrameRotation(2,1)<<
" "<<
577 globalToAlignFrameRotation(2,2));
580 const int nAlignPar = alignPars->
size();
584 for(
int i(0);
i<nAlignPar+3; ++
i) derivatives[
i].setZero();
588 for (; iatsos != alignTrack->
lastAtsos(); ++iatsos) {
595 int nResDim = alignTSOS->
nResDim();
602 std::vector<Amg::VectorX> * atsosDerivs =
nullptr;
603 std::vector<Amg::VectorX> * atsosDerVtx =
nullptr;
605 atsosDerivs =
new std::vector<Amg::VectorX>(nResDim,
Amg::VectorX(nAlignPar));
606 atsosDerVtx =
new std::vector<Amg::VectorX>(nResDim,
Amg::VectorX(3));
607 ATH_MSG_DEBUG(
"nResDim = "<<nResDim<<
" vector size is "<<atsosDerivs->size());
608 ATH_MSG_DEBUG(
"nAlignPar = "<<nAlignPar<<
" CLHEP::HepVector size is "<<atsosDerivs->at(0).rows());
615 if (!mtp || !(mtp->covariance()) )
continue;
621 localToGlobalRotation(0,1) <<
" " <<
622 localToGlobalRotation(0,2));
624 localToGlobalRotation(1,1) <<
" " <<
625 localToGlobalRotation(1,2));
627 localToGlobalRotation(2,1) <<
" " <<
628 localToGlobalRotation(2,2));
630 if(
double alphastrip=alignTSOS->
alphaStrip()) {
631 ATH_MSG_DEBUG(
"applying fanout rotation : " << alphastrip );
635 localToGlobalRotation(0,1) <<
" " <<
636 localToGlobalRotation(0,2));
638 localToGlobalRotation(1,1) <<
" " <<
639 localToGlobalRotation(1,2));
641 localToGlobalRotation(2,1) <<
" " <<
642 localToGlobalRotation(2,2));
679 ATH_MSG_DEBUG(
"trackdir " << trackdir[0] <<
" " << trackdir[1] <<
" " << trackdir[2]);
683 double cotphi_x = trackdir.x() / trackdir.z();
688 cotphi_x = trackdir.y() / trackdir.z();
691 double Rxx = R(0,0) - cotphi_x * R(0,2);
692 double Ryx = R(1,0) - cotphi_x * R(1,2);
693 double Rzx = R(2,0) - cotphi_x * R(2,2);
694 ATH_MSG_DEBUG(
"Rxx/Ryx/Rzx: " << Rxx <<
"/" << Ryx <<
"/" << Rzx);
730 const double z0z0 = 366.5*366.5;
737 Amg::Vector3D RxGlob=-1.0 * (globalToAlignFrameRotation.inverse() * RxLoc);
739 for (
int ipar=0; ipar<nAlignPar; ipar++) {
740 const AlignPar * alignPar = (*alignPars)[ipar];
744 derivatives[ipar][imeas] = projR[paramType];
746 (*atsosDerivs)[0][ipar] = projR[paramType];
749 for (
int ipar=0; ipar<3; ipar++) {
750 derivatives[nAlignPar+ipar][imeas] = RxGlob[ipar];
754 for (
int i=0;
i<nAlignPar+3;
i++)
755 ATH_MSG_DEBUG(
"derivatives["<<
i<<
"]["<<imeas<<
"]="<<derivatives[
i][imeas]);
762 double cotphi_y = trackdir.y() / trackdir.z() ;
763 double Rxy = R(0,1) - cotphi_y * R(0,2) ;
764 double Ryy = R(1,1) - cotphi_y * R(1,2) ;
765 double Rzy = R(2,1) - cotphi_y * R(2,2) ;
766 ATH_MSG_DEBUG(
"Rxy/Ryy/Rzy: " << Rxy <<
"/" << Ryy <<
"/" << Rzy);
783 Amg::Vector3D RyGlob=-1.0 * (globalToAlignFrameRotation.inverse() * RyLoc);
785 for (
int ipar=0; ipar<nAlignPar; ipar++) {
786 const AlignPar * alignPar = (*alignPars)[ipar];
788 ATH_MSG_DEBUG(
"2nd dim, ipar="<<ipar<<
", paramType="<<paramType);
790 derivatives[ipar][imeas] = projR[paramType];
792 (*atsosDerivs)[1][ipar] = projR[paramType];
796 for (
int ipar=0; ipar<3; ipar++) {
797 derivatives[nAlignPar+ipar][imeas] = RyGlob[ipar];
801 for (
int i=0;
i<nAlignPar+3;
i++)
802 ATH_MSG_DEBUG(
"2nd dim: derivatives["<<
i<<
"]["<<imeas<<
"]="<<derivatives[
i][imeas]);