Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:56:02 

0001 //   COCOA class implementation file
0002 //Id:  LightRay.cc
0003 //CAT: Fit
0004 //
0005 //   History: v1.0
0006 //   Pedro Arce
0007 
0008 #include "Alignment/CocoaModel/interface/LightRay.h"
0009 #include "Alignment/CocoaModel/interface/ALIPlane.h"
0010 #include "Alignment/CocoaUtilities/interface/ALIUtils.h"
0011 #include "Alignment/CocoaModel/interface/OpticalObject.h"
0012 #include <cstdlib>
0013 #include <cmath>  // include floating-point std::abs functions
0014 
0015 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0016 //@@ construct a default LightRay
0017 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0018 LightRay::LightRay() {
0019   _point = CLHEP::Hep3Vector(0., 0., 0.);
0020   _direction = CLHEP::Hep3Vector(0., 0., 1.);
0021 }
0022 
0023 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0024 //@@ Set the point and direction to that of the laser or source
0025 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0026 void LightRay::startLightRay(OpticalObject* opto) {
0027   if (ALIUtils::debug >= 3)
0028     std::cout << std::endl << "LR: CREATE LIGHTRAY " << opto->name() << " OptO type is " << opto->type() << std::endl;
0029 
0030   //---------- Get Z axis of opto
0031   CLHEP::Hep3Vector ZAxis(0., 0., 1.);
0032   const CLHEP::HepRotation& rmt = opto->rmGlob();
0033   ZAxis = rmt * ZAxis;
0034 
0035   //---------- By convention, direction of LightRay = opto_ZAxis
0036   setDirection(ZAxis);
0037   setPoint(opto->centreGlob());
0038 
0039   if (ALIUtils::debug >= 3) {
0040     dumpData(" LightRay at creation ");
0041   }
0042   if (ALIUtils::debug >= 5) {
0043     ALIUtils::dumprm(rmt, "laser Rotation matrix");
0044   }
0045 }
0046 
0047 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0048 LightRay::LightRay(OpticalObject* opto1, OpticalObject* opto2) {
0049   if (ALIUtils::debug >= 7)
0050     std::cout << std::endl << "LR:CREATE LIGHTRAY FROM SOURCE" << opto2->name() << std::endl;
0051 
0052   CLHEP::Hep3Vector _ZAxis(0., 0., 1.);
0053   //-  LightRay* linetmp;
0054   //-linetmp = new LightRay;
0055   //---------- set direction and point
0056   setDirection(opto2->centreGlob() - opto1->centreGlob());
0057   setPoint(opto1->centreGlob());
0058 
0059   if (ALIUtils::debug >= 9)
0060     std::cout << "OPT" << opto1 << opto1->name() << std::endl;
0061   //-  std::cout << "centre glob" << &p1->aff()->centre_glob() << std::endl;
0062   if (ALIUtils::debug >= 9) {
0063     dumpData(" ");
0064   }
0065 }
0066 
0067 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0068 LightRay::LightRay(CLHEP::Hep3Vector& vec1, CLHEP::Hep3Vector& vec2) {
0069   CLHEP::Hep3Vector dir = vec2 - vec1;
0070   dir *= 1. / dir.mag();
0071   setDirection(dir);
0072   setPoint(vec1);
0073   if (ALIUtils::debug >= 9) {
0074     dumpData(" ");
0075   }
0076 }
0077 
0078 //@@ intersect: Intersect a LightRay with a plane and change thePoint to the intersection point
0079 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0080 void LightRay::intersect(const OpticalObject& opto) {
0081   if (ALIUtils::debug >= 3)
0082     std::cout << "% LR INTERSECT WITH OPTO" << std::endl;
0083   CLHEP::Hep3Vector ZAxis(0., 0., 1.);
0084   const CLHEP::HepRotation& rmt = opto.rmGlob();
0085   ZAxis = rmt * ZAxis;
0086   ALIPlane optoPlane(opto.centreGlob(), ZAxis);
0087   intersect(optoPlane);
0088 }
0089 
0090 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0091 //@@ intersect: Intersect a LightRay with a plane and change thePoint to the intersection point
0092 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0093 void LightRay::intersect(const ALIPlane& plane) {
0094   if (ALIUtils::debug >= 4)
0095     std::cout << "% LR INTERSECT WITH PLANE" << std::endl;
0096   if (ALIUtils::debug >= 4) {
0097     ALIUtils::dump3v(plane.point(), "plane point");
0098     ALIUtils::dump3v(plane.normal(), "plane normal");
0099     //-    dumpData(" ");
0100   }
0101 
0102   //---------- Check that they intersect
0103   if (std::abs(plane.normal() * direction()) < 1.E-10) {
0104     std::cerr << " !!!! INTERSECTION NOT POSSIBLE: LightRay is perpendicular to plane " << std::endl;
0105     std::cerr << " plane.normal()*direction() = " << plane.normal() * direction() << std::endl;
0106     ALIUtils::dump3v(direction(), "LightRay direction ");
0107     ALIUtils::dump3v(plane.normal(), "plane normal ");
0108     exit(1);
0109   }
0110 
0111   //---------- Get intersection point between LightRay and plane
0112   CLHEP::Hep3Vector vtemp = plane.point() - _point;
0113   if (ALIUtils::debug >= 5)
0114     ALIUtils::dump3v(vtemp, "n_r = point  - point_plane");
0115   ALIdouble dtemp = _direction * plane.normal();
0116   if (ALIUtils::debug >= 5)
0117     std::cout << " lightray* plate normal" << dtemp << std::endl;
0118   if (dtemp != 0.) {
0119     dtemp = (vtemp * plane.normal()) / dtemp;
0120     if (ALIUtils::debug >= 5)
0121       std::cout << " n_r*plate normal" << dtemp << std::endl;
0122   } else {
0123     std::cerr << "!!! LightRay: Intersect With Plane: plane and light ray parallel: no intersection" << std::endl;
0124   }
0125   vtemp = _direction * dtemp;
0126   if (ALIUtils::debug >= 5)
0127     ALIUtils::dump3v(vtemp, "n_r scaled");
0128   CLHEP::Hep3Vector inters = vtemp + _point;
0129   if (ALIUtils::debug >= 3)
0130     ALIUtils::dump3v(inters, "INTERSECTION point ");
0131 
0132   _point = inters;
0133 }
0134 
0135 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0136 //@@ Intersect the LightRay with a plane and then change the direction from reflection on this plane
0137 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0138 void LightRay::reflect(const ALIPlane& plane) {
0139   intersect(plane);
0140   if (ALIUtils::debug >= 4)
0141     std::cout << "% LR: REFLECT IN PLANE " << std::endl;
0142   ALIdouble cosang = -(plane.normal() * _direction) / plane.normal().mag() / _direction.mag();
0143   if (ALIUtils::debug >= 5) {
0144     std::cout << " cosang = " << cosang << std::endl;
0145     ALIUtils::dump3v(plane.normal(), " plane normal");
0146   }
0147   _direction += plane.normal() * 2 * cosang;
0148   if (ALIUtils::debug >= 5)
0149     dumpData("LightRay after reflection: ");
0150 }
0151 
0152 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0153 //@@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
0154 //@@ 3D deviation: it actually deviates in the plane of the plate normal and lightray, in the other plane there is no deviation
0155 //@@ Refract inside this plane
0156 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0157 void LightRay::refract(const ALIPlane& plate, const ALIdouble refra_ind1, const ALIdouble refra_ind2) {
0158   if (ALIUtils::debug >= 5) {
0159     std::cout << "% LR REFRACT: "
0160               << "refra_ind1 = " << refra_ind1 << " refra_ind2 = " << refra_ind2 << std::endl;
0161     std::cout << "@ First intersect with plate plane " << std::endl;
0162   }
0163 
0164   intersect(plate);
0165 
0166   //---------- First plane: formed by plate normal and lightray, but get two ortonormal std::vectors in this plane (one of it plate normal)
0167   CLHEP::Hep3Vector Axis1 = plate.normal().cross(direction());
0168   //----- Check lightray is not parallel to plate normal
0169   if (Axis1.mag() < 1.E-6) {
0170     if (ALIUtils::debug >= 3) {
0171       std::cout << " light ray normal to plane, no refraction " << std::endl;
0172     }
0173     if (ALIUtils::debug >= 2) {
0174       dumpData("LightRay after refraction: ");
0175     }
0176 
0177     return;
0178   }
0179 
0180   if (ALIUtils::debug >= 5) {
0181     ALIUtils::dump3v(Axis1, " axis 1 temp ");
0182   }
0183   Axis1 = Axis1.cross(plate.normal());
0184   Axis1 *= 1. / Axis1.mag();
0185   //----- Project lightray on this plane
0186   if (ALIUtils::debug >= 4) {
0187     ALIUtils::dump3v(plate.normal(), " plate normal ");
0188     ALIUtils::dump3v(Axis1, " axis 1 ");
0189   }
0190 
0191   //----- Angle between LightRay and plate_normal before traversing
0192   ALIdouble cosang = -(plate.normal() * direction()) / plate.normal().mag() / direction().mag();
0193   ALIdouble sinang = sqrt(1. - cosang * cosang);
0194 
0195   //----- Angle between LightRay projection and plate normal after traversing (refracted)
0196   ALIdouble sinangp = sinang * refra_ind1 / refra_ind2;
0197   if (std::abs(sinangp) > 1.) {
0198     std::cerr << " !!!EXITING LightRay::refract: incidence ray on plane too close to face, refraction will not allow "
0199                  "entering "
0200               << std::endl;
0201     ALIUtils::dump3v(plate.normal(), " plate normal ");
0202     ALIUtils::dump3v(direction(), " light ray direction ");
0203     std::cout << " refraction index first medium " << refra_ind1 << " refraction index second medium " << refra_ind2
0204               << std::endl;
0205     exit(1);
0206   }
0207 
0208   if (ALIUtils::debug >= 4) {
0209     std::cout << "LightRay refract on plane 1: sin(ang) before = " << sinang << " sinang after " << sinangp
0210               << std::endl;
0211   }
0212   ALIdouble cosangp = sqrt(1. - sinangp * sinangp);
0213   //----- Change Lightray direction in this plane
0214   //--- Get sign of projections in plate normal and axis1
0215   ALIdouble signN = direction() * plate.normal();
0216   signN /= std::abs(signN);
0217   ALIdouble sign1 = direction() * Axis1;
0218   sign1 /= std::abs(sign1);
0219   if (ALIUtils::debug >= 4) {
0220     dumpData("LightRay refract: direction before plate");
0221     std::cout << " sign projection on plate normal " << signN << " sign projection on Axis1 " << sign1 << std::endl;
0222   }
0223   setDirection(signN * cosangp * plate.normal() + sign1 * sinangp * Axis1);
0224   //-  std::cout << " " << signN  << " " << cosangp  << " " << plate.normal() << " " << sign1  << " " << sinangp   << " " << Axis1 << std::endl;
0225 
0226   if (ALIUtils::debug >= 3) {
0227     dumpData("LightRay refract: direction after plate");
0228   }
0229 }
0230 
0231 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0232 //@@ shift and deviates around X, Y and Z of opto
0233 //@@
0234 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0235 void LightRay::shiftAndDeviateWhileTraversing(const OpticalObject* opto, char behav) {
0236   ALIstring ename("devi X");
0237   ename[4] = behav;
0238   ename[5] = 'X';
0239   ALIdouble deviX = opto->findExtraEntryValue(ename);
0240   ename[5] = 'Y';
0241   ALIdouble deviY = opto->findExtraEntryValue(ename);
0242   ename[5] = 'Z';
0243   ALIdouble deviZ = opto->findExtraEntryValue(ename);
0244 
0245   ename = "shift X";
0246   ename[5] = behav;
0247   ename[6] = 'X';
0248   ALIdouble shiftX = opto->findExtraEntryValue(ename);
0249   ename[6] = 'Y';
0250   ALIdouble shiftY = opto->findExtraEntryValue(ename);
0251   ename[6] = 'Z';
0252   ALIdouble shiftZ = opto->findExtraEntryValue(ename);
0253 
0254   if (ALIUtils::debug >= 3) {
0255     //-    std::cout << " shift X " << shiftX << " shiftY " << shiftY << " shiftZ " << shiftZ << std::endl;
0256     //-    std::cout << " deviX " << deviX << " deviY " << deviY << " deviZ " << deviZ << std::endl;
0257     std::cout << " shift X " << shiftX << " shift Y " << shiftY << std::endl;
0258     std::cout << " devi X " << deviX << " devi Y " << deviY << std::endl;
0259   }
0260 
0261   shiftAndDeviateWhileTraversing(opto, shiftX, shiftY, shiftZ, deviX, deviY, deviZ);
0262   //  shiftAndDeviateWhileTraversing( shiftX, shiftY, deviX, deviY );
0263 }
0264 
0265 void LightRay::shiftAndDeviateWhileTraversing(const OpticalObject* opto,
0266                                               ALIdouble shiftX,
0267                                               ALIdouble shiftY,
0268                                               ALIdouble shiftZ,
0269                                               ALIdouble deviX,
0270                                               ALIdouble deviY,
0271                                               ALIdouble deviZ) {
0272   //----- Get local opto X, Y and Z axis
0273   CLHEP::Hep3Vector XAxis(1., 0., 0.);
0274   CLHEP::Hep3Vector YAxis(0., 1., 0.);
0275   CLHEP::Hep3Vector ZAxis(0., 0., 1.);
0276   const CLHEP::HepRotation& rmt = opto->rmGlob();
0277   XAxis = rmt * XAxis;
0278   YAxis = rmt * YAxis;
0279   ZAxis = rmt * ZAxis;
0280 
0281   if (ALIUtils::debug >= 5) {
0282     ALIUtils::dump3v(XAxis, "X axis of opto");
0283     ALIUtils::dump3v(YAxis, "Y axis of opto");
0284     ALIUtils::dump3v(ZAxis, "Z axis of opto");
0285   }
0286 
0287   //---------- Shift
0288   CLHEP::Hep3Vector pointold = _point;
0289   _point += shiftX * XAxis;
0290   _point += shiftY * YAxis;
0291   _point += shiftZ * ZAxis;
0292   if (_point != pointold && ALIUtils::debug >= 3) {
0293     ALIUtils::dump3v(_point - pointold, "CHANGE point");
0294   }
0295 
0296   //---------- Deviate
0297   CLHEP::Hep3Vector direcold = _direction;
0298   if (ALIUtils::debug >= 5) {
0299     ALIUtils::dump3v(XAxis, "XAxis");
0300     ALIUtils::dump3v(YAxis, "YAxis");
0301     ALIUtils::dump3v(ZAxis, "ZAxis");
0302     ALIUtils::dump3v(_direction, "LightRay direction");
0303   }
0304 
0305   _direction.rotate(deviX, XAxis);
0306   if (_direction != direcold && ALIUtils::debug >= 3) {
0307     std::cout << " deviX " << deviX << std::endl;
0308     ALIUtils::dump3v(_direction - direcold, "CHANGE direction");
0309   }
0310   _direction.rotate(deviY, YAxis);
0311   if (_direction != direcold && ALIUtils::debug >= 3) {
0312     std::cout << " deviY " << deviY << std::endl;
0313     ALIUtils::dump3v(_direction - direcold, "CHANGE direction");
0314   }
0315   _direction.rotate(deviZ, ZAxis);
0316   if (_direction != direcold && ALIUtils::debug >= 3) {
0317     std::cout << " deviZ " << deviZ << std::endl;
0318     ALIUtils::dump3v(_direction - direcold, "CHANGE direction");
0319   }
0320 
0321   if (_direction != direcold && ALIUtils::debug >= 3) {
0322     ALIUtils::dump3v(_direction - direcold, "CHANGE direction");
0323   }
0324 }
0325 
0326 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0327 //@@ shift and deviates around two directions perpendicular to LightRay
0328 //@@
0329 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0330 /*void LightRay::shiftAndDeviateWhileTraversing( ALIdouble shiftAxis1, ALIdouble shiftAxis2, ALIdouble deviAxis1, ALIdouble deviAxis2 )
0331 {
0332   if(ALIUtils::debug >= 3) {
0333     std::cout << " shift Axis 1 " << shiftAxis1 << " shift Axis 2 " << shiftAxis2 << std::endl;
0334     std::cout << " devi Axis 1 " << deviAxis1 << " devi Axis 2 " << deviAxis2 << std::endl;
0335   }
0336 
0337   //----- Get two directions perpendicular to LightRay
0338   //-- First can be (y,-x,0), unless direciton is (0,0,1), or close
0339   CLHEP::Hep3Vector PerpAxis1;
0340   if(std::abs(std::abs(_direction.z())-1.) > 0.1) {
0341     if (ALIUtils::debug >= 99) ALIUtils::dump3v( _direction, "_direction1");
0342     PerpAxis1 = CLHEP::Hep3Vector(_direction.y(),-_direction.x(),0.);
0343   } else {
0344     if (ALIUtils::debug >= 99) ALIUtils::dump3v( _direction, "_direction2");
0345     PerpAxis1 = CLHEP::Hep3Vector(_direction.z(),0.,-_direction.y());
0346   }
0347   if (ALIUtils::debug >= 4) ALIUtils::dump3v( PerpAxis1, "1st perpendicular direction of DDet");
0348 
0349   CLHEP::Hep3Vector PerpAxis2 = _direction.cross(PerpAxis1);
0350   if (ALIUtils::debug >= 4) ALIUtils::dump3v( PerpAxis2, "2nd perpendicular direction of DDet");
0351 
0352   //---------- Shift
0353   CLHEP::Hep3Vector pointold = _point;
0354   _point += shiftAxis1*PerpAxis1;
0355   _point += shiftAxis2*PerpAxis2;
0356   if(_point != pointold && ALIUtils::debug >= 3 ) {
0357     ALIUtils::dump3v( _point-pointold, "CHANGE point");
0358   }
0359 
0360   //---------- Deviate
0361   CLHEP::Hep3Vector direcold = _direction;
0362   _direction.rotate(deviAxis1, PerpAxis1);
0363   _direction.rotate(deviAxis2, PerpAxis2);
0364   if(_direction != direcold && ALIUtils::debug >= 3) {
0365     ALIUtils::dump3v( _direction-direcold, "CHANGE direction");
0366   }
0367 
0368 }
0369 */
0370 
0371 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0372 //@@
0373 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0374 void LightRay::dumpData(const ALIstring& str) const {
0375   std::cout << str << std::endl;
0376   ALIUtils::dump3v(_point, "$$ LightRay point: ");
0377   ALIUtils::dump3v(_direction, "$$ LightRay direction: ");
0378   /*
0379   CLHEP::Hep3Vector dirn = _direction;
0380   dirn.rotateZ( -23.72876*3.1415926/180.);
0381   ALIUtils::dump3v( dirn, "$$ LightRay direction: ");
0382   */
0383 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!