CgemGeoAlign.cxx

Go to the documentation of this file.

#include "CgemGeomSvc/CgemGeoAlign.h"
 
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/StatusCode.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/PropertyMgr.h"
 
#include "TMath.h"
 
#include <iostream>
#include <fstream>
#include <iomanip>
 
using namespace std;
 
CgemGeoAlign::CgemGeoAlign(){
     m_r[0] = 87.915;
     m_r[1] = 130.415;
     m_r[2] = 172.915;
}
 
void CgemGeoAlign::initAlignPar(string alignFile){
     ifstream fin(alignFile.c_str());
     char strtmp[200];
     fin.getline(strtmp,1000);
     string str[3];
     for (int layer=0; layer<3; layer++){
      fin >> str[layer] >> m_dx[layer] >> m_dy[layer] >> m_dz[layer]
          >> m_rx[layer] >> m_ry[layer] >> m_rz[layer];
 
      m_dxOrig[layer] = m_dx[layer];
      m_dyOrig[layer] = m_dy[layer];
      m_dzOrig[layer] = m_dz[layer];
      m_rxOrig[layer] = m_rx[layer];
      m_ryOrig[layer] = m_ry[layer];
      m_rzOrig[layer] = m_rz[layer];
     }
     fin.close();
}
 
void CgemGeoAlign::resetAlignPar(){
     for (int layer=0; layer<3; layer++){
      m_dx[layer] = m_dxOrig[layer];
      m_dy[layer] = m_dyOrig[layer];
      m_dz[layer] = m_dzOrig[layer];
      m_rx[layer] = m_rxOrig[layer];
      m_ry[layer] = m_ryOrig[layer];
      m_rz[layer] = m_rzOrig[layer];
     }
}
 
StraightLine CgemGeoAlign::StraightLineConversion(int layer, StraightLine lineOriginal){
     // 1.  Get 2 points on straight line:
     HepPoint3D p1 = lineOriginal.xAtR(m_r[layer], -1);
     HepPoint3D p2 = lineOriginal.xAtR(m_r[layer], 1);
 
     double dist = p1.distance2(p2);
     if(fabs(dist) < 0.001){
      double s1 = lineOriginal.sAtR(m_r[layer], -1);
      double s2 = lineOriginal.sAtR(m_r[layer], 1);
 
      p1 = lineOriginal.x(s1);
      p2 = lineOriginal.x(s2);
     }
 
     // 2.  Transformation
     // cout << "alignment par: " << setw(15) << m_dx[layer] << setw(15) << m_dy[layer]
     //       << setw(15) << m_dz[layer] << setw(15) << m_rz[layer] << endl;
 
     // cout << "2 points before transform: "
     //       << setw(15) << p1.x() << setw(15) << p1.y() << setw(15) << p1.z()
     //       << setw(15) << p2.x() << setw(15) << p2.y() << setw(15) << p2.z() << endl;
 
     HepPoint3D newP1 = point_transform(layer, p1);
     HepPoint3D newP2 = point_transform(layer, p2);
 
     // cout << "2 points after transform:  "
     //       << setw(15) << newP1.x() << setw(15) << newP1.y() << setw(15) << newP1.z()
     //       << setw(15) << newP2.x() << setw(15) << newP2.y() << setw(15) << newP2.z() << endl;
 
     // 3.  New line parameters
     StraightLine newLine(newP1, newP2);
 
     // cout << "use StraightLine" << endl;
     // cout << "track before transform: " << setw(15) << lineOriginal.dr() << setw(15) << lineOriginal.phi0()
     //       << setw(15) << lineOriginal.dz() << setw(15) << lineOriginal.tanl() << endl;
     // cout << "track after transform:  " << setw(15) << newLine.dr() << setw(15) << newLine.phi0()
     //       << setw(15) << newLine.dz() << setw(15) << newLine.tanl() << endl << endl;
 
     return newLine;
}
 
void CgemGeoAlign::StraightLineConversion_v1(int layer, double lineOriginal[], double lineConverted[]){
 
     // m_dx[0] = 1;
     // m_dy[0] = 1;
     // m_dz[0] = 1;
     // m_rz[0] = 0; //TMath::Pi()/2.;
 
     // 0. Get original straight line parameters:
          
     double drho = lineOriginal[0];
     double phi0 = lineOriginal[1];
     double dz   = lineOriginal[2];
     double tgl  = lineOriginal[3];
 
     // 1.  Get 2 points on straight line:
 
     //--- Line on x-y plane: y = ax + b
     //--- Line on s-z plane: z = -1.*tgl*s + dz;
     //                       s = sqrt((x-x0)^2+(y-y0)^2) 
 
     int flg_parallel_x = 0;
     int flg_parallel_y = 0;
 
     if(phi0 == TMath::Pi()/2. || phi0 == -1.*TMath::Pi()/2.)     flg_parallel_x = 1;
     if(phi0 == 0 || phi0 == TMath::Pi())                         flg_parallel_y = 1;
 
     double a, b, x0, x1, x2, y0, y1, y2, s1, s2, z1, z2;
     x0 = drho*cos(phi0);
     y0 = drho*sin(phi0);
 
     if(flg_parallel_x == 0 && flg_parallel_y == 0){
      a  = tan(TMath::Pi()/2.+phi0);
      b  = drho / cos(phi0);
 
      x1 = 0;
      y1 = drho / cos(phi0);
      s1 = sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
      z1 = -1.*tgl*s1 + dz;
 
      x2 = -1.*b/a;
      y2 = 0;
      s2 = sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0));
      z2 = -1.*tgl*s2 + dz;
      if(x1==x2 || y1==y2 || z1==z2){
           x2 = (1.-b)/a;
           y2 = 1;
           s2 = sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0));
           z2 = -1.*tgl*s2 + dz;
      }
     }
     else if(flg_parallel_x == 1){
      x1 = drho; y1 = 0; s1 = sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0)); z1 = -1.*tgl*s1 + dz;
      x2 = drho; y2 = 1; s2 = sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0)); z2 = -1.*tgl*s2 + dz;
     }
     else if(flg_parallel_y == 1){
      x1 = 0; y1 = drho; s1 = sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0)); z1 = -1.*tgl*s1 + dz;
      x2 = 1; y2 = drho; s2 = sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0)); z2 = -1.*tgl*s2 + dz;
     }
 
     // 2.  Transformation
 
     double shift_x    = m_dx[layer]; 
     double shift_y    = m_dy[layer]; 
     double shift_z    = m_dz[layer];
     double rotation_z = m_rz[layer];
 
     // cout << "alignment par: " << setw(15) << shift_x << setw(15) << shift_y
     //       << setw(15) << shift_z << setw(15) << rotation_z << endl;
 
     // cout << "2 points before transform: "
     //       << setw(15) << x1 << setw(15) << y1 << setw(15) << z1
     //       << setw(15) << x2 << setw(15) << y2 << setw(15) << z2 << endl;
 
     point_transform(x1,  y1,  z1,  shift_x,  shift_y,  shift_z,  rotation_z);
     point_transform(x2,  y2,  z2,  shift_x,  shift_y,  shift_z,  rotation_z);
 
     // cout << "2 points after transform:  "
     //       << setw(15) << x1 << setw(15) << y1 << setw(15) << z1
     //       << setw(15) << x2 << setw(15) << y2 << setw(15) << z2 << endl;
 
     // 3.  New line parameters
 
     //--- Line on x-y plane: y = ax + b
     //          y1   = ax1 + b
     //          y2   = ax2 + b
     //          a    = (y1-y2)/(x1-x2) = -tan(TMath::Pi()/2.-phi0) = tan(phi0-TMath::Pi()/2.);
     //          b    = y1-a*x1 =  drho / cos(phi0);
     //          phi0 = atan(a)+TMath::Pi()/2.;
     //          drho = b*cos(phi0);
     //--- Line on s-z plane: z = -1.*tgl*s + dz;
     //                       s = sqrt((x-x0)^2+(y-y0)^2)
     //          z1   = -1.*tgl*sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))+dz 
     //          z2   = -1.*tgl*sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0))+dz
     //          tgl  = (z1-z2)/(-1.*(sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))-sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0))));
     //          dz   = z1 + tgl*sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
 
     double new_a;   
     double new_b;   
     double new_phi0;
     double new_drho;
     double new_tgl; 
     double new_dz;  
 
     if(x1!=x2 && y1!=y2){
      new_a    = (y1-y2)/(x1-x2);  
      new_b    = y1-new_a*x1;      
      new_phi0 = atan(new_a)+TMath::Pi()/2.;
      new_drho = new_b*cos(phi0);
      new_tgl  = (z1-z2)/(-1.*(sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))-sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0))));
      new_dz   = z1 + new_tgl*sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
     }
     else if(x1==x2){
      new_phi0 = 0;
      new_drho = x1;
      new_tgl  = (z1-z2)/(-1.*(sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))-sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0))));
      new_dz   = z1 + new_tgl*sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
     }
     else if(y1==y2){
      new_phi0 = TMath::Pi()/2.;
      new_drho = y1;
      new_tgl  = (z1-z2)/(-1.*(sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0))-sqrt((x2-x0)*(x2-x0)+(y2-y0)*(y2-y0))));
      new_dz   = z1 + new_tgl*sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
     }
 
     lineConverted[0] = new_drho;
     lineConverted[1] = new_phi0;
     lineConverted[2] = new_dz;  
     lineConverted[3] = new_tgl; 
 
     // cout << "track before transform: " << setw(15) << lineOriginal[0] << setw(15) << lineOriginal[1]
     //       << setw(15) << lineOriginal[2] << setw(15) << lineOriginal[3] << endl;
     // cout << "track after transform:  " << setw(15) << lineConverted[0] << setw(15) << lineConverted[1]
     //       << setw(15) << lineConverted[2] << setw(15) << lineConverted[3] << endl << endl;
 
}
 
void CgemGeoAlign::HelixConversion(int layer, double helixOriginal[], double helixConverted[]){
}
 
HepPoint3D CgemGeoAlign::point_transform(int layer, HepPoint3D pos){
     double x_shift = pos.x() - m_dx[layer];
     double y_shift = pos.y() - m_dy[layer];
     double z_shift = pos.z() - m_dz[layer];
 
     double x =     x_shift*cos(m_rz[layer]) + y_shift*sin(m_rz[layer]);
     double y = -1.*x_shift*sin(m_rz[layer]) + y_shift*cos(m_rz[layer]);
     double z =     z_shift;
 
     return HepPoint3D(x,y,z);
}
 
void CgemGeoAlign::point_transform(double& x, double& y, double& z, double shift_x,
                   double shift_y, double shift_z, double rotation_z){
     // double x_origin = x;
     // double y_origin = y;
     // double z_origin = z;
 
     // double x_rotate =     x_origin*cos(rotation_z) + y_origin*sin(rotation_z);
     // double y_rotate = -1.*x_origin*sin(rotation_z) + y_origin*cos(rotation_z);
     // double z_rotate =     z_origin;
 
     // double x_shift = x_rotate + shift_x;
     // double y_shift = y_rotate + shift_y;
     // double z_shift = z_rotate + shift_z;
 
     // x = x_shift;
     // y = y_shift;
     // z = z_shift;
 
     double x_shift = x - shift_x;
     double y_shift = y - shift_y;
     double z_shift = z - shift_z;
 
     x =     x_shift*cos(rotation_z) + y_shift*sin(rotation_z);
     y = -1.*x_shift*sin(rotation_z) + y_shift*cos(rotation_z);
     z =     z_shift;
}
 
HepPoint3D CgemGeoAlign::point_invTransform(int layer, HepPoint3D pos){
     double xp = pos.x();
     double yp = pos.y();
     double zp = pos.z();
 
     double rz = -1.*m_rz[layer];
     double xrot =     xp*cos(rz) + yp*sin(rz);
     double yrot = -1.*xp*sin(rz) + yp*cos(rz);
 
     double x = xrot + m_dx[layer];
     double y = yrot + m_dy[layer];
     double z = zp + m_dz[layer];
 
     // cout << "before invTransform: " << setw(15) << xp<< setw(15) << yp << setw(15) << zp << endl;
     // cout << " after invTransform: " << setw(15) << xp<< setw(15) << yp << setw(15) << zp << endl;
 
     return HepPoint3D(x, y, z);
}