LightRay.cc

Go to the documentation of this file.

//   COCOA class implementation file
//Id:  LightRay.cc
//CAT: Fit
//
//   History: v1.0
//   Pedro Arce

#include "Alignment/CocoaModel/interface/LightRay.h"
#include "Alignment/CocoaModel/interface/ALIPlane.h"
#include "Alignment/CocoaUtilities/interface/ALIUtils.h"
#include "Alignment/CocoaModel/interface/OpticalObject.h"
#include <cstdlib>
#include <cmath>        // include floating-point std::abs functions


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//@@ construct a default LightRay
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
LightRay::LightRay( )
{
  _point = CLHEP::Hep3Vector(0.,0.,0.);
  _direction = CLHEP::Hep3Vector(0.,0.,1.);
}

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//@@ Set the point and direction to that of the laser or source
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void LightRay::startLightRay( OpticalObject* opto )
{
  if(ALIUtils::debug >= 3) std::cout << std::endl << "LR: CREATE LIGHTRAY " << opto->name() << " OptO type is " << opto->type() << std::endl;

  //---------- Get Z axis of opto
  CLHEP::Hep3Vector ZAxis(0.,0.,1.);
  const CLHEP::HepRotation& rmt = opto->rmGlob();
  ZAxis = rmt * ZAxis;

  //---------- By convention, direction of LightRay = opto_ZAxis
  setDirection( ZAxis );
  setPoint( opto->centreGlob() );

 if (ALIUtils::debug >= 3) {
    dumpData(" LightRay at creation ");
  }
  if (ALIUtils::debug >= 5) {
    ALIUtils::dumprm( rmt, "laser Rotation matrix");
  }

}


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
LightRay::LightRay( OpticalObject* opto1, OpticalObject* opto2)
{
  if(ALIUtils::debug >= 7) std::cout << std::endl << "LR:CREATE LIGHTRAY FROM SOURCE" << opto2->name() << std::endl;

  CLHEP::Hep3Vector _ZAxis(0.,0.,1.);
  //-  LightRay* linetmp;
  //-linetmp = new LightRay;
  //---------- set direction and point
  setDirection( opto2->centreGlob() - opto1->centreGlob() );
  setPoint( opto1->centreGlob() );

  if(ALIUtils::debug >= 9) std::cout << "OPT" << opto1 << opto1->name() << std::endl;
  //-  std::cout << "centre glob" << &p1->aff()->centre_glob() << std::endl;
  if (ALIUtils::debug >= 9) {
    dumpData(" ");
  }

}


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
LightRay::LightRay( CLHEP::Hep3Vector& vec1, CLHEP::Hep3Vector& vec2 )
{
  CLHEP::Hep3Vector dir = vec2 - vec1;
  dir *= 1./dir.mag();
  setDirection( dir );
  setPoint( vec1 );
  if (ALIUtils::debug >= 9) {
    dumpData(" ");
  }
}


//@@ intersect: Intersect a LightRay with a plane and change thePoint to the intersection point
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void LightRay::intersect( const OpticalObject& opto )
{
  if(ALIUtils::debug >= 3) std::cout << "% LR INTERSECT WITH OPTO" << std::endl;
  CLHEP::Hep3Vector ZAxis(0.,0.,1.);
  const CLHEP::HepRotation& rmt = opto.rmGlob();
  ZAxis = rmt*ZAxis;
  ALIPlane optoPlane(opto.centreGlob(), ZAxis);
  intersect( optoPlane );
}


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//@@ intersect: Intersect a LightRay with a plane and change thePoint to the intersection point
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void LightRay::intersect( const ALIPlane& plane )
{
  if(ALIUtils::debug >= 4) std::cout << "% LR INTERSECT WITH PLANE" << std::endl;
  if(ALIUtils::debug >= 4) {
    ALIUtils::dump3v( plane.point(), "plane point");
    ALIUtils::dump3v( plane.normal(), "plane normal");
    //-    dumpData(" ");
  }

  //---------- Check that they intersect
  if( std::abs( plane.normal()*direction() ) < 1.E-10 ) {
    std::cerr << " !!!! INTERSECTION NOT POSSIBLE: LightRay is perpendicular to plane " << std::endl;
    std::cerr << " plane.normal()*direction() = " << plane.normal()*direction() << std::endl;
    ALIUtils::dump3v( direction(), "LightRay direction ");
    ALIUtils::dump3v( plane.normal(), "plane normal ");
    exit(1);
  }

  //---------- Get intersection point between LightRay and plane
  CLHEP::Hep3Vector vtemp = plane.point() - _point;
  if(ALIUtils::debug >= 5) ALIUtils::dump3v( vtemp, "n_r = point  - point_plane");
  ALIdouble dtemp = _direction * plane.normal();
  if(ALIUtils::debug >= 5) std::cout << " lightray* plate normal" << dtemp << std::endl;
  if ( dtemp != 0. ) {
      dtemp = (vtemp * plane.normal()) / dtemp;
      if(ALIUtils::debug >= 5)  std::cout << " n_r*plate normal" << dtemp << std::endl;
  } else {
      std::cerr << "!!! LightRay: Intersect With Plane: plane and light ray parallel: no intersection" << std::endl;
  }
  vtemp = _direction * dtemp;
  if(ALIUtils::debug >= 5) ALIUtils::dump3v( vtemp, "n_r scaled");
  CLHEP::Hep3Vector inters = vtemp + _point;
  if(ALIUtils::debug >= 3) ALIUtils::dump3v( inters, "INTERSECTION point ");

  _point = inters;
}


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//@@ Intersect the LightRay with a plane and then change the direction from reflection on this plane
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void LightRay::reflect( const ALIPlane& plane)
{
  intersect(plane);
  if( ALIUtils::debug >= 4 ) std::cout << "% LR: REFLECT IN PLANE " << std::endl;
  ALIdouble cosang = -(plane.normal() * _direction) /
           plane.normal().mag() / _direction.mag();
  if(ALIUtils::debug >= 5) {
    std::cout << " cosang = " << cosang << std::endl;
    ALIUtils::dump3v( plane.normal(), " plane normal");
  }
  _direction += plane.normal()*2*cosang;
  if( ALIUtils::debug >= 5 ) dumpData("LightRay after reflection: ");

}


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//@@Deviate a LightRay because of refraction when it passes from a medium of refraction index  refra_ind1 to a medium of refraction index  refra_ind2
//@@ 3D deviation: it actually deviates in the plane of the plate normal and lightray, in the other plane there is no deviation
//@@ Refract inside this plane
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void LightRay::refract( const ALIPlane& plate, const ALIdouble refra_ind1, const ALIdouble refra_ind2)
{
  if(ALIUtils::debug >= 5) {
    std::cout << "% LR REFRACT: " <<  "refra_ind1 = " << refra_ind1 << " refra_ind2 = " << refra_ind2 <<std::endl;
    std::cout << "@ First intersect with plate plane " << std::endl;
  }

  intersect( plate );

  //---------- First plane: formed by plate normal and lightray, but get two ortonormal std::vectors in this plane (one of it plate normal)
  CLHEP::Hep3Vector Axis1 = plate.normal().cross( direction() );
  //----- Check lightray is not parallel to plate normal
  if( Axis1.mag() < 1.E-6 ) {
    if(ALIUtils::debug >= 3) {
      std::cout << " light ray normal to plane, no refraction " << std::endl;
    }
    if (ALIUtils::debug >= 2) {
      dumpData("LightRay after refraction: ");
    }

    return;
  }

  if(ALIUtils::debug >= 5) {
    ALIUtils::dump3v( Axis1,  " axis 1 temp ");
  }
  Axis1 = Axis1.cross( plate.normal() );
  Axis1 *= 1./Axis1.mag();
  //----- Project lightray on this plane
  if(ALIUtils::debug >= 4) {
    ALIUtils::dump3v( plate.normal(),  " plate normal ");
    ALIUtils::dump3v( Axis1,  " axis 1 ");
  }

  //----- Angle between LightRay and plate_normal before traversing
  ALIdouble cosang = -(plate.normal() * direction()) /
                     plate.normal().mag() / direction().mag();
  ALIdouble sinang = sqrt( 1. - cosang*cosang );

  //----- Angle between LightRay projection and plate normal after traversing (refracted)
  ALIdouble sinangp = sinang * refra_ind1 / refra_ind2;
  if( std::abs(sinangp) > 1. ) {
    std::cerr << " !!!EXITING LightRay::refract: incidence ray on plane too close to face, refraction will not allow entering " << std::endl;
    ALIUtils::dump3v( plate.normal(), " plate normal ");
    ALIUtils::dump3v( direction(), " light ray direction ");
    std::cout << " refraction index first medium " << refra_ind1  << " refraction index second medium " << refra_ind2 << std::endl;
    exit(1);
  }

  if(ALIUtils::debug >= 4) {
    std::cout << "LightRay refract on plane 1: sin(ang) before = " << sinang
     << " sinang after " << sinangp << std::endl;
  }
  ALIdouble cosangp = sqrt( 1. - sinangp*sinangp );
  //----- Change Lightray direction in this plane
  //--- Get sign of projections in plate normal and axis1
  ALIdouble signN = direction()*plate.normal();
  signN /= std::abs(signN);
  ALIdouble sign1 = direction()*Axis1;
  sign1 /= std::abs(sign1);
  if(ALIUtils::debug >= 4) {
    dumpData("LightRay refract: direction before plate");
    std::cout << " sign projection on plate normal " << signN << " sign projection on Axis1 " << sign1 << std::endl;
  }
  setDirection( signN * cosangp * plate.normal() + sign1 * sinangp * Axis1);
  //-  std::cout << " " << signN  << " " << cosangp  << " " << plate.normal() << " " << sign1  << " " << sinangp   << " " << Axis1 << std::endl;

  if(ALIUtils::debug >= 3) {
    dumpData("LightRay refract: direction after plate");
  }

}


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//@@ shift and deviates around X, Y and Z of opto
//@@
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void LightRay::shiftAndDeviateWhileTraversing( const OpticalObject* opto, char behav )
{
  ALIstring ename("devi X");
  ename[4] = behav;
  ename[5] = 'X';
  ALIdouble deviX = opto->findExtraEntryValue(ename);
  ename[5] = 'Y';
  ALIdouble deviY = opto->findExtraEntryValue(ename);
  ename[5] = 'Z';
  ALIdouble deviZ = opto->findExtraEntryValue(ename);

  ename = "shift X";
  ename[5] = behav;
  ename[6] = 'X';
  ALIdouble shiftX = opto->findExtraEntryValue(ename);
  ename[6] = 'Y';
  ALIdouble shiftY = opto->findExtraEntryValue(ename);
  ename[6] = 'Z';
  ALIdouble shiftZ = opto->findExtraEntryValue(ename);

  if(ALIUtils::debug >= 3) {
    //-    std::cout << " shift X " << shiftX << " shiftY " << shiftY << " shiftZ " << shiftZ << std::endl;
    //-    std::cout << " deviX " << deviX << " deviY " << deviY << " deviZ " << deviZ << std::endl;
    std::cout << " shift X " << shiftX << " shift Y " << shiftY << std::endl;
    std::cout << " devi X " << deviX << " devi Y " << deviY << std::endl;
  }

  shiftAndDeviateWhileTraversing( opto, shiftX, shiftY, shiftZ, deviX, deviY, deviZ );
  //  shiftAndDeviateWhileTraversing( shiftX, shiftY, deviX, deviY );
}


void LightRay::shiftAndDeviateWhileTraversing( const OpticalObject* opto, ALIdouble shiftX, ALIdouble shiftY, ALIdouble shiftZ, ALIdouble deviX, ALIdouble deviY, ALIdouble deviZ )
{

  //----- Get local opto X, Y and Z axis
  CLHEP::Hep3Vector XAxis(1.,0.,0.);
  CLHEP::Hep3Vector YAxis(0.,1.,0.);
  CLHEP::Hep3Vector ZAxis(0.,0.,1.);
  const CLHEP::HepRotation& rmt = opto->rmGlob();
  XAxis = rmt*XAxis;
  YAxis = rmt*YAxis;
  ZAxis = rmt*ZAxis;

  if (ALIUtils::debug >= 5) {
    ALIUtils::dump3v( XAxis, "X axis of opto");
    ALIUtils::dump3v( YAxis, "Y axis of opto");
    ALIUtils::dump3v( ZAxis, "Z axis of opto");
  }

  //---------- Shift
  CLHEP::Hep3Vector pointold = _point;
  _point += shiftX*XAxis;
  _point += shiftY*YAxis;
  _point += shiftZ*ZAxis;
  if(_point != pointold && ALIUtils::debug >= 3 ) {
    ALIUtils::dump3v( _point-pointold, "CHANGE point");
  }

  //---------- Deviate
  CLHEP::Hep3Vector direcold = _direction;
  if( ALIUtils::debug >= 5) {
    ALIUtils::dump3v( XAxis, "XAxis");
    ALIUtils::dump3v( YAxis, "YAxis");
    ALIUtils::dump3v( ZAxis, "ZAxis");
    ALIUtils::dump3v( _direction, "LightRay direction");
  }

  _direction.rotate(deviX, XAxis);
  if(_direction != direcold && ALIUtils::debug >= 3) {
    std::cout << " deviX " << deviX << std::endl;
    ALIUtils::dump3v( _direction-direcold, "CHANGE direction");
  }
  _direction.rotate(deviY, YAxis);
  if(_direction != direcold && ALIUtils::debug >= 3) {
    std::cout << " deviY " << deviY << std::endl;
    ALIUtils::dump3v( _direction-direcold, "CHANGE direction");
  }
  _direction.rotate(deviZ, ZAxis);
  if(_direction != direcold && ALIUtils::debug >= 3) {
    std::cout << " deviZ " << deviZ << std::endl;
    ALIUtils::dump3v( _direction-direcold, "CHANGE direction");
  }

  if(_direction != direcold && ALIUtils::debug >= 3) {
    ALIUtils::dump3v( _direction-direcold, "CHANGE direction");
  }

}



//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//@@ shift and deviates around two directions perpendicular to LightRay
//@@
 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
/*void LightRay::shiftAndDeviateWhileTraversing( ALIdouble shiftAxis1, ALIdouble shiftAxis2, ALIdouble deviAxis1, ALIdouble deviAxis2 )
{
  if(ALIUtils::debug >= 3) {
    std::cout << " shift Axis 1 " << shiftAxis1 << " shift Axis 2 " << shiftAxis2 << std::endl;
    std::cout << " devi Axis 1 " << deviAxis1 << " devi Axis 2 " << deviAxis2 << std::endl;
  }

  //----- Get two directions perpendicular to LightRay
  //-- First can be (y,-x,0), unless direciton is (0,0,1), or close
  CLHEP::Hep3Vector PerpAxis1;
  if(std::abs(std::abs(_direction.z())-1.) > 0.1) {
    if (ALIUtils::debug >= 99) ALIUtils::dump3v( _direction, "_direction1");
    PerpAxis1 = CLHEP::Hep3Vector(_direction.y(),-_direction.x(),0.);
  } else {
    if (ALIUtils::debug >= 99) ALIUtils::dump3v( _direction, "_direction2");
    PerpAxis1 = CLHEP::Hep3Vector(_direction.z(),0.,-_direction.y());
  }
  if (ALIUtils::debug >= 4) ALIUtils::dump3v( PerpAxis1, "1st perpendicular direction of DDet");

  CLHEP::Hep3Vector PerpAxis2 = _direction.cross(PerpAxis1);
  if (ALIUtils::debug >= 4) ALIUtils::dump3v( PerpAxis2, "2nd perpendicular direction of DDet");

  //---------- Shift
  CLHEP::Hep3Vector pointold = _point;
  _point += shiftAxis1*PerpAxis1;
  _point += shiftAxis2*PerpAxis2;
  if(_point != pointold && ALIUtils::debug >= 3 ) {
    ALIUtils::dump3v( _point-pointold, "CHANGE point");
  }

  //---------- Deviate
  CLHEP::Hep3Vector direcold = _direction;
  _direction.rotate(deviAxis1, PerpAxis1);
  _direction.rotate(deviAxis2, PerpAxis2);
  if(_direction != direcold && ALIUtils::debug >= 3) {
    ALIUtils::dump3v( _direction-direcold, "CHANGE direction");
  }

}
*/


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//@@
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void LightRay::dumpData(const ALIstring& str) const
{
  std::cout << str << std::endl;
  ALIUtils::dump3v( _point, "$$ LightRay point: ");
  ALIUtils::dump3v( _direction, "$$ LightRay direction: ");
  /*
  CLHEP::Hep3Vector dirn = _direction;
  dirn.rotateZ( -23.72876*3.1415926/180.);
  ALIUtils::dump3v( dirn, "$$ LightRay direction: ");
  */
}