---

PointForAlignment.cc

Go to the documentation of this file.

#include "Alignment/MuonStandaloneAlgorithm/interface/PointForAlignment.h"
#include "TMath.h"

PointForAlignment::PointForAlignment(){}

void PointForAlignment::setPointForAlignment(long rawId_, int subDet_M, int station_M, GlobalPoint globalPosition_M, GlobalVector globalDirection_M, LocalPoint localPosition_M, LocalError localPositionError_M, LocalVector localDirection_M, LocalError localDirectionError_M, GlobalPoint globalPosition_P, GlobalVector globalDirection_P, LocalPoint localPosition_P, LocalVector localDirection_P,TMatrixD theDerivatives_, TMatrixD CovMatrix_m) {
  
  //The rawId of the detector
  theRawId = rawId_;
  subDet = subDet_M;
  station = station_M;
  
  //Position and orientation vectors
  theGlobalPosition = globalPosition_M;
  theGlobalDirection = globalDirection_M;
  theLocalPosition = localPosition_M;
  theLocalPositionError = localPositionError_M;
  theLocalDirection = localDirection_M;
  theLocalDirectionError = localDirectionError_M;
  //Position and orientation predicted vectors
  thePredictedGlobalPosition = globalPosition_P;
  thePredictedGlobalDirection = globalDirection_P;
  thePredictedLocalPosition = localPosition_P;
  thePredictedLocalDirection = localDirection_P;
  
  //Jacobian
  theDerivatives.ResizeTo(theDerivatives_.GetNrows(), theDerivatives_.GetNcols());  
  theDerivatives = theDerivatives_;
  
  //Covariance Matrix
  CovMatrixInit.ResizeTo(CovMatrix_m.GetNrows(), CovMatrix_m.GetNcols());
  CovMatrixInit = CovMatrix_m;

  calculateErrors();
  calculateResiduals();
  calculateAlignmentCoeff();

  if(TMath::Abs(theResiduals(0,0)) < 5.0 && TMath::Abs(theResiduals(1,0)) < 8.0 && TMath::Abs(theResiduals(2,0)) < 1.0 && TMath::Abs(theResiduals(3,0) < 1.0)) {
    validPoint = true;
  } else {
    validPoint = false;
  } 
  
  if(TMath::Abs(theDerivatives(1,0)/theDerivatives(0,0)) > 0.1 ||
     TMath::Abs(theDerivatives(1,1)/theDerivatives(0,1)) > 0.1 ||
     TMath::Abs(theDerivatives(1,3)/theDerivatives(0,3)) > 0.1 ||
     TMath::Abs(theDerivatives(0,2)/theDerivatives(1,2)) > 0.1 ||
     TMath::Abs(theDerivatives(0,4)/theDerivatives(1,4)) > 0.1) validPoint = false; 
  }

PointForAlignment::~PointForAlignment(){}

long PointForAlignment::rawId(){return theRawId;}

GlobalPoint PointForAlignment::globalPosition(){return theGlobalPosition;}

GlobalVector PointForAlignment::globalDirection(){return theGlobalDirection;}

LocalPoint PointForAlignment::localPosition(){return theLocalPosition;}

LocalVector PointForAlignment::localDirection(){return theLocalDirection;}

GlobalPoint PointForAlignment::predictedGlobalPosition(){return thePredictedGlobalPosition;}

GlobalVector PointForAlignment::predictedGlobalDirection(){return thePredictedGlobalDirection;}

LocalPoint PointForAlignment::predictedLocalPosition(){return thePredictedLocalPosition;}

LocalVector PointForAlignment::predictedLocalDirection(){return thePredictedLocalDirection;}

TMatrixD PointForAlignment::derivatives(){return theDerivatives;}
 
TMatrixD PointForAlignment::residual() {return theResiduals;}
 
TMatrixD PointForAlignment::errors() {return theErrors;}

TMatrixD PointForAlignment::alignmentMatrix() {return theAlignmentMatrix;}

double PointForAlignment::residualRPhi(){return theResidualRPhi;}

double PointForAlignment::residualZ(){return theResidualZ;}

double PointForAlignment::residualPhi(){return theResidualPhi;}

double PointForAlignment::residualTheta(){return theResidualTheta;}


void PointForAlignment::calculateResiduals() {
  
  //Calculation of residuals
  if(station != 4 && subDet == 1) {
    theResidualRPhi = localPosition().x() - predictedLocalPosition().x();
    theResidualZ = localPosition().y()-predictedLocalPosition().y();
    theResidualPhi = TMath::ATan2(localDirection().z(),localDirection().x())-TMath::ATan2(predictedLocalDirection().z(),predictedLocalDirection().x());
    theResidualTheta = TMath::ATan2(localDirection().z(), localDirection().y())-TMath::ATan2(predictedLocalDirection().z(), predictedLocalDirection().y());
  } else {
    theResidualRPhi = localPosition().x() - predictedLocalPosition().x();
    theResidualZ = 0;
    theResidualPhi = TMath::ATan2(localDirection().z(),localDirection().x())-TMath::ATan2(predictedLocalDirection().z(),predictedLocalDirection().x());
    theResidualTheta = 0;
  }
  
  theResiduals.ResizeTo(NDOFCoor,1);
  theResiduals(0,0) = residualRPhi();
  theResiduals(1,0) = residualZ();
  theResiduals(2,0) = residualPhi();
  theResiduals(3,0) = residualTheta();

}

void PointForAlignment::calculateErrors() {
 
  theErrors.ResizeTo(NDOFCoor, NDOFCoor);
  CovMatrix.ResizeTo(NDOFCoor, NDOFCoor);
  TMatrixD unit(NDOFCoor, NDOFCoor);
  unit.UnitMatrix();
  TMatrixD comparison(NDOFCoor, NDOFCoor);
  double crosscheck;
  double det;
  if(station != 4 && subDet == 1) {
    CovMatrix = CovMatrixInit;
    theErrors = CovMatrix;
    theErrors.Invert(&det);
    comparison = unit - CovMatrix*theErrors;
    crosscheck = comparison.E2Norm();
  } else {
    CovMatrix.ResizeTo(NDOFCoor/2, NDOFCoor/2);
    TMatrixD checkMatrix(NDOFCoor/2, NDOFCoor/2);
    CovMatrix(0,0) = CovMatrixInit(0,0);
    CovMatrix(0,1) = CovMatrixInit(0,2);
    CovMatrix(1,0) = CovMatrixInit(2,0);
    CovMatrix(1,1) = CovMatrixInit(2,2);
    checkMatrix = CovMatrix;
    CovMatrix.Invert(&det); 
    comparison.ResizeTo(NDOFCoor/2, NDOFCoor/2);
    unit.ResizeTo(NDOFCoor/2, NDOFCoor/2);
    unit.UnitMatrix();
    theErrors(0,0) = CovMatrix(0,0);
    theErrors(0,2) = CovMatrix(0,1);
    theErrors(2,0) = CovMatrix(1,0);
    theErrors(2,2) = CovMatrix(1,1);
    comparison = unit - CovMatrix*checkMatrix;
    crosscheck = comparison.E2Norm();
  }
  if(TMath::Abs(det) < 1.e-30 || crosscheck > 1.0e-15) {
    validPoint = false;
    std::cout << "<PointForAlignment> Not valid point" << std::endl;
  }
}

void PointForAlignment::calculateAlignmentCoeff() {

  double tanphi = localDirection().z()/localDirection().x();
  double phi = TMath::ATan2(localDirection().z(),localDirection().x());
  double cosphi = TMath::Cos(phi);
  double sinphi = TMath::Sin(phi);
  double tantheta = localDirection().z()/localDirection().y();
  double theta = TMath::ATan2(localDirection().z(),localDirection().y());
  double costheta = TMath::Cos(theta);
  double sintheta = TMath::Sin(theta);
  theAlignmentMatrix.ResizeTo(NDOFCoor, NDOFAlign);
  if(station != 4 && subDet == 1) {
    if(tanphi == 0.0 || tantheta == 0.0) {
      validPoint = false;
      return;
    }
    //X
    theAlignmentMatrix(0,0) = 1.0; 
    theAlignmentMatrix(0,1) = 0.0;
    theAlignmentMatrix(0,2) = -1.0/tanphi; 
    theAlignmentMatrix(0,3) = -localPosition().y();
    theAlignmentMatrix(0,4) = localPosition().x()/tanphi;
    theAlignmentMatrix(0,5) = -1.0/tanphi*localPosition().y();
    //Z
    theAlignmentMatrix(1,0) = 0.0;
    theAlignmentMatrix(1,1) = 1.0;
    theAlignmentMatrix(1,2) = -1.0/tantheta;
    theAlignmentMatrix(1,3) = localPosition().x();
    theAlignmentMatrix(1,4) = localPosition().x()/tantheta;
    theAlignmentMatrix(1,5) = -localPosition().y()/tantheta;
    //Phi
    theAlignmentMatrix(2,0) = 0.0; theAlignmentMatrix(2,1) = 0.0;
    theAlignmentMatrix(2,2) = 0.0; theAlignmentMatrix(2,3) = 0.0;
    theAlignmentMatrix(2,4) = cosphi*cosphi; theAlignmentMatrix(2,5) = cosphi*sinphi/tantheta;
    //theta
    theAlignmentMatrix(3,0) = 0.0; theAlignmentMatrix(3,1) = 0.0;
    theAlignmentMatrix(3,2) = 0.0; theAlignmentMatrix(3,3) = 0.0;
    theAlignmentMatrix(3,4) = costheta*sintheta/tanphi; theAlignmentMatrix(3,5) = costheta*costheta;
  } else {
    if(tanphi == 0.0) {
      validPoint = false;
      return;
    }
    //X
    theAlignmentMatrix(0,0) = 1.0; 
    theAlignmentMatrix(0,1) = 0.0;
    theAlignmentMatrix(0,2) = -1.0/tanphi; 
    theAlignmentMatrix(0,3) = 0.0; 
    theAlignmentMatrix(0,4) = localPosition().x()/tanphi;
    theAlignmentMatrix(0,5) = 0.0;
    //Z
    theAlignmentMatrix(1,0) = 0.0;
    theAlignmentMatrix(1,1) = 0.0;
    theAlignmentMatrix(1,2) = 0.0;
    theAlignmentMatrix(1,3) = 0.0;
    theAlignmentMatrix(1,4) = 0.0;
    theAlignmentMatrix(1,5) = 0.0;
    //Phi
    theAlignmentMatrix(2,0) = 0.0; theAlignmentMatrix(2,1) = 0.0;
    theAlignmentMatrix(2,2) = 0.0; theAlignmentMatrix(2,3) = 0.0;
    theAlignmentMatrix(2,4) = cosphi*cosphi; theAlignmentMatrix(2,5) = 0.0;
    //theta
    theAlignmentMatrix(3,0) = 0.0; theAlignmentMatrix(3,1) = 0.0;
    theAlignmentMatrix(3,2) = 0.0; theAlignmentMatrix(3,3) = 0.0;
    theAlignmentMatrix(3,4) = 0.0; theAlignmentMatrix(3,5) = 0.0;
  }
}


bool PointForAlignment::pointIsValid() {return validPoint;}

---