Reference Manual

Go to the documentation of this file.
00001 //   COCOA class implementation file
00002 //Id:  LightRay.cc
00003 //CAT: Fit
00004 //
00005 //   History: v1.0 
00006 //   Pedro Arce
00007 
00008 #include "Alignment/CocoaModel/interface/LightRay.h"
00009 #include "Alignment/CocoaModel/interface/ALIPlane.h" 
00010 #include "Alignment/CocoaUtilities/interface/ALIUtils.h"
00011 #include "Alignment/CocoaModel/interface/OpticalObject.h"
00012 #include <cstdlib>
00013 
00014 
00015 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00016 //@@ construct a default LightRay
00017 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00018 LightRay::LightRay( )
00019 {
00020   _point = CLHEP::Hep3Vector(0.,0.,0.);
00021   _direction = CLHEP::Hep3Vector(0.,0.,1.);
00022 }
00023 
00024 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00025 //@@ Set the point and direction to that of the laser or source
00026 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00027 void LightRay::startLightRay( OpticalObject* opto )
00028 {
00029   if(ALIUtils::debug >= 3) std::cout << std::endl << "LR: CREATE LIGHTRAY " << opto->name() << " OptO type is " << opto->type() << std::endl;
00030 
00031   //---------- Get Z axis of opto
00032   CLHEP::Hep3Vector ZAxis(0.,0.,1.);
00033   CLHEP::HepRotation rmt = opto->rmGlob();
00034   ZAxis = rmt * ZAxis;
00035 
00036   //---------- By convention, direction of LightRay = opto_ZAxis
00037   setDirection( ZAxis );
00038   setPoint( opto->centreGlob() );
00039 
00040  if (ALIUtils::debug >= 3) {
00041     dumpData(" LightRay at creation ");
00042   }
00043   if (ALIUtils::debug >= 5) {
00044     ALIUtils::dumprm( rmt, "laser Rotation matrix");
00045   }
00046 
00047 } 
00048 
00049 
00050 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00051 LightRay::LightRay( OpticalObject* opto1, OpticalObject* opto2)
00052 {
00053   if(ALIUtils::debug >= 7) std::cout << std::endl << "LR:CREATE LIGHTRAY FROM SOURCE" << opto2->name() << std::endl;
00054 
00055   CLHEP::Hep3Vector _ZAxis(0.,0.,1.);
00056   //-  LightRay* linetmp;
00057   //-linetmp = new LightRay; 
00058   //---------- set direction and point
00059   setDirection( opto2->centreGlob() - opto1->centreGlob() );
00060   setPoint( opto1->centreGlob() );
00061 
00062   if(ALIUtils::debug >= 9) std::cout << "OPT" << opto1 << opto1->name() << std::endl;
00063   //-  std::cout << "centre glob" << &p1->aff()->centre_glob() << std::endl;
00064   if (ALIUtils::debug >= 9) {
00065     dumpData(" ");
00066   }
00067 
00068 }
00069 
00070 
00071 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00072 LightRay::LightRay( CLHEP::Hep3Vector& vec1, CLHEP::Hep3Vector& vec2 )
00073 {
00074   CLHEP::Hep3Vector dir = vec2 - vec1;
00075   dir *= 1./dir.mag();
00076   setDirection( dir );
00077   setPoint( vec1 );
00078   if (ALIUtils::debug >= 9) {
00079     dumpData(" ");
00080   }
00081 }
00082 
00083 
00084 //@@ intersect: Intersect a LightRay with a plane and change thePoint to the intersection point
00085 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00086 void LightRay::intersect( const OpticalObject& opto )
00087 {
00088   if(ALIUtils::debug >= 3) std::cout << "% LR INTERSECT WITH OPTO" << std::endl;
00089   CLHEP::Hep3Vector ZAxis(0.,0.,1.);
00090   CLHEP::HepRotation rmt = opto.rmGlob(); 
00091   ZAxis = rmt*ZAxis;
00092   ALIPlane optoPlane(opto.centreGlob(), ZAxis);
00093   intersect( optoPlane );
00094 }
00095 
00096 
00097 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00098 //@@ intersect: Intersect a LightRay with a plane and change thePoint to the intersection point
00099 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00100 void LightRay::intersect( const ALIPlane& plane )
00101 {
00102   if(ALIUtils::debug >= 4) std::cout << "% LR INTERSECT WITH PLANE" << std::endl;
00103   if(ALIUtils::debug >= 4) {
00104     ALIUtils::dump3v( plane.point(), "plane point"); 
00105     ALIUtils::dump3v( plane.normal(), "plane normal"); 
00106     //-    dumpData(" ");
00107   }
00108   
00109   //---------- Check that they intersect
00110   if( fabs( plane.normal()*direction() ) < 1.E-10 ) {
00111     std::cerr << " !!!! INTERSECTION NOT POSSIBLE: LightRay is perpendicular to plane " << std::endl;
00112     std::cerr << " plane.normal()*direction() = " << plane.normal()*direction() << std::endl;
00113     ALIUtils::dump3v( direction(), "LightRay direction ");
00114     ALIUtils::dump3v( plane.normal(), "plane normal ");
00115     exit(1);
00116   }
00117 
00118   //---------- Get intersection point between LightRay and plane
00119   CLHEP::Hep3Vector vtemp = plane.point() - _point;
00120   if(ALIUtils::debug >= 5) ALIUtils::dump3v( vtemp, "n_r = point  - point_plane"); 
00121   ALIdouble dtemp = _direction * plane.normal();
00122   if(ALIUtils::debug >= 5) std::cout << " lightray* plate normal" << dtemp << std::endl;
00123   if ( dtemp != 0. ) {
00124       dtemp = (vtemp * plane.normal()) / dtemp;
00125       if(ALIUtils::debug >= 5)  std::cout << " n_r*plate normal" << dtemp << std::endl;
00126   } else {
00127       std::cerr << "!!! LightRay: Intersect With Plane: plane and light ray parallel: no intersection" << std::endl;     
00128   } 
00129   vtemp = _direction * dtemp;
00130   if(ALIUtils::debug >= 5) ALIUtils::dump3v( vtemp, "n_r scaled"); 
00131   CLHEP::Hep3Vector inters = vtemp + _point;
00132   if(ALIUtils::debug >= 3) ALIUtils::dump3v( inters, "INTERSECTION point "); 
00133 
00134   _point = inters;
00135 }
00136 
00137 
00138 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00139 //@@ Intersect the LightRay with a plane and then change the direction from reflection on this plane
00140 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00141 void LightRay::reflect( const ALIPlane& plane)
00142 {
00143   intersect(plane);
00144   if( ALIUtils::debug >= 4 ) std::cout << "% LR: REFLECT IN PLANE " << std::endl;
00145   ALIdouble cosang = -(plane.normal() * _direction) / 
00146            plane.normal().mag() / _direction.mag();
00147   if(ALIUtils::debug >= 5) { 
00148     std::cout << " cosang = " << cosang << std::endl;
00149     ALIUtils::dump3v( plane.normal(), " plane normal");
00150   }
00151   _direction += plane.normal()*2*cosang;
00152   if( ALIUtils::debug >= 5 ) dumpData("LightRay after reflection: ");
00153 
00154 }
00155 
00156 
00157 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00158 //@@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
00159 //@@ 3D deviation: it actually deviates in the plane of the plate normal and lightray, in the other plane there is no deviation
00160 //@@ Refract inside this plane
00161 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00162 void LightRay::refract( const ALIPlane plate, const ALIdouble refra_ind1, const ALIdouble refra_ind2)
00163 {
00164   if(ALIUtils::debug >= 5) {
00165     std::cout << "% LR REFRACT: " <<  "refra_ind1 = " << refra_ind1 << " refra_ind2 = " << refra_ind2 <<std::endl;
00166     std::cout << "@ First intersect with plate plane " << std::endl;
00167   }
00168 
00169   intersect( plate );
00170 
00171   //---------- First plane: formed by plate normal and lightray, but get two ortonormal std::vectors in this plane (one of it plate normal)
00172   CLHEP::Hep3Vector Axis1 = plate.normal().cross( direction() );
00173   //----- Check lightray is not parallel to plate normal
00174   if( Axis1.mag() < 1.E-6 ) {
00175     if(ALIUtils::debug >= 3) {
00176       std::cout << " light ray normal to plane, no refraction " << std::endl;
00177     }
00178     if (ALIUtils::debug >= 2) {
00179       dumpData("LightRay after refraction: "); 
00180     }
00181 
00182     return;
00183   }
00184 
00185   if(ALIUtils::debug >= 5) {
00186     ALIUtils::dump3v( Axis1,  " axis 1 temp ");
00187   }
00188   Axis1 = Axis1.cross( plate.normal() );
00189   Axis1 *= 1./Axis1.mag();
00190   //----- Project lightray on this plane 
00191   if(ALIUtils::debug >= 4) {
00192     ALIUtils::dump3v( plate.normal(),  " plate normal ");
00193     ALIUtils::dump3v( Axis1,  " axis 1 ");
00194   }
00195 
00196   //----- Angle between LightRay and plate_normal before traversing
00197   ALIdouble cosang = -(plate.normal() * direction()) / 
00198                      plate.normal().mag() / direction().mag();
00199   ALIdouble sinang = sqrt( 1. - cosang*cosang );
00200   
00201   //----- Angle between LightRay projection and plate normal after traversing (refracted)
00202   ALIdouble sinangp = sinang * refra_ind1 / refra_ind2; 
00203   if( fabs(sinangp) > 1. ) {
00204     std::cerr << " !!!EXITING LightRay::refract: incidence ray on plane too close to face, refraction will not allow entering " << std::endl;
00205     ALIUtils::dump3v( plate.normal(), " plate normal ");
00206     ALIUtils::dump3v( direction(), " light ray direction ");
00207     std::cout << " refraction index first medium " << refra_ind1  << " refraction index second medium " << refra_ind2 << std::endl; 
00208     exit(1);
00209   }
00210 
00211   if(ALIUtils::debug >= 4) {
00212     std::cout << "LightRay refract on plane 1: sin(ang) before = " << sinang 
00213          << " sinang after " << sinangp << std::endl;
00214   }
00215   ALIdouble cosangp = sqrt( 1. - sinangp*sinangp );
00216   //----- Change Lightray direction in this plane 
00217   //--- Get sign of projections in plate normal and axis1
00218   ALIdouble signN = direction()*plate.normal();
00219   signN /= fabs(signN);
00220   ALIdouble sign1 = direction()*Axis1;
00221   sign1 /= fabs(sign1);
00222   if(ALIUtils::debug >= 4) {
00223     dumpData("LightRay refract: direction before plate");
00224     std::cout << " sign projection on plate normal " << signN << " sign projection on Axis1 " << sign1 << std::endl;
00225   }
00226   setDirection( signN * cosangp * plate.normal() + sign1 * sinangp * Axis1);
00227   //-  std::cout << " " << signN  << " " << cosangp  << " " << plate.normal() << " " << sign1  << " " << sinangp   << " " << Axis1 << std::endl;
00228   
00229   if(ALIUtils::debug >= 3) {
00230     dumpData("LightRay refract: direction after plate");
00231   }
00232 
00233 }
00234 
00235 
00236 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00237 //@@ shift and deviates around X, Y and Z of opto
00238 //@@  
00239 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00240 void LightRay::shiftAndDeviateWhileTraversing( const OpticalObject* opto, char behav )
00241 {
00242   ALIstring ename("devi X");
00243   ename[4] = behav;
00244   ename[5] = 'X';
00245   ALIdouble deviX = opto->findExtraEntryValue(ename);
00246   ename[5] = 'Y';
00247   ALIdouble deviY = opto->findExtraEntryValue(ename);
00248   ename[5] = 'Z';
00249   ALIdouble deviZ = opto->findExtraEntryValue(ename);
00250   
00251   ename = "shift X";
00252   ename[5] = behav;
00253   ename[6] = 'X';
00254   ALIdouble shiftX = opto->findExtraEntryValue(ename);
00255   ename[6] = 'Y';
00256   ALIdouble shiftY = opto->findExtraEntryValue(ename);
00257   ename[6] = 'Z';
00258   ALIdouble shiftZ = opto->findExtraEntryValue(ename);
00259 
00260   if(ALIUtils::debug >= 3) {
00261     //-    std::cout << " shift X " << shiftX << " shiftY " << shiftY << " shiftZ " << shiftZ << std::endl;
00262     //-    std::cout << " deviX " << deviX << " deviY " << deviY << " deviZ " << deviZ << std::endl;
00263     std::cout << " shift X " << shiftX << " shift Y " << shiftY << std::endl;
00264     std::cout << " devi X " << deviX << " devi Y " << deviY << std::endl;
00265   }
00266 
00267   shiftAndDeviateWhileTraversing( opto, shiftX, shiftY, shiftZ, deviX, deviY, deviZ );
00268   //  shiftAndDeviateWhileTraversing( shiftX, shiftY, deviX, deviY );
00269 } 
00270 
00271  
00272 void LightRay::shiftAndDeviateWhileTraversing( const OpticalObject* opto, ALIdouble shiftX, ALIdouble shiftY, ALIdouble shiftZ, ALIdouble deviX, ALIdouble deviY, ALIdouble deviZ )
00273 {
00274 
00275   //----- Get local opto X, Y and Z axis
00276   CLHEP::Hep3Vector XAxis(1.,0.,0.);
00277   CLHEP::Hep3Vector YAxis(0.,1.,0.);
00278   CLHEP::Hep3Vector ZAxis(0.,0.,1.);
00279   CLHEP::HepRotation rmt = opto->rmGlob();
00280   XAxis = rmt*XAxis;
00281   YAxis = rmt*YAxis;
00282   ZAxis = rmt*ZAxis;
00283   
00284   if (ALIUtils::debug >= 5) {
00285     ALIUtils::dump3v( XAxis, "X axis of opto");
00286     ALIUtils::dump3v( YAxis, "Y axis of opto");
00287     ALIUtils::dump3v( ZAxis, "Z axis of opto");
00288   }
00289 
00290   //---------- Shift
00291   CLHEP::Hep3Vector pointold = _point;
00292   _point += shiftX*XAxis;
00293   _point += shiftY*YAxis;
00294   _point += shiftZ*ZAxis;
00295   if(_point != pointold && ALIUtils::debug >= 3 ) {
00296     ALIUtils::dump3v( _point-pointold, "CHANGE point");
00297   }
00298 
00299   //---------- Deviate  
00300   CLHEP::Hep3Vector direcold = _direction;
00301   if( ALIUtils::debug >= 5) { 
00302     ALIUtils::dump3v( XAxis, "XAxis");
00303     ALIUtils::dump3v( YAxis, "YAxis");
00304     ALIUtils::dump3v( ZAxis, "ZAxis");
00305     ALIUtils::dump3v( _direction, "LightRay direction");
00306   }
00307   
00308   _direction.rotate(deviX, XAxis);
00309   if(_direction != direcold && ALIUtils::debug >= 3) {
00310     std::cout << " deviX " << deviX << std::endl;
00311     ALIUtils::dump3v( _direction-direcold, "CHANGE direction");
00312   }
00313   _direction.rotate(deviY, YAxis);
00314   if(_direction != direcold && ALIUtils::debug >= 3) {
00315     std::cout << " deviY " << deviY << std::endl;
00316     ALIUtils::dump3v( _direction-direcold, "CHANGE direction");
00317   }
00318   _direction.rotate(deviZ, ZAxis);
00319   if(_direction != direcold && ALIUtils::debug >= 3) {
00320     std::cout << " deviZ " << deviZ << std::endl;
00321     ALIUtils::dump3v( _direction-direcold, "CHANGE direction");
00322   }
00323 
00324   if(_direction != direcold && ALIUtils::debug >= 3) {
00325     ALIUtils::dump3v( _direction-direcold, "CHANGE direction");
00326   }
00327 
00328 }
00329 
00330 
00331 
00332 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00333 //@@ shift and deviates around two directions perpendicular to LightRay
00334 //@@  
00335  //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00336 /*void LightRay::shiftAndDeviateWhileTraversing( ALIdouble shiftAxis1, ALIdouble shiftAxis2, ALIdouble deviAxis1, ALIdouble deviAxis2 )
00337 {
00338   if(ALIUtils::debug >= 3) {
00339     std::cout << " shift Axis 1 " << shiftAxis1 << " shift Axis 2 " << shiftAxis2 << std::endl;
00340     std::cout << " devi Axis 1 " << deviAxis1 << " devi Axis 2 " << deviAxis2 << std::endl;
00341   }
00342 
00343   //----- Get two directions perpendicular to LightRay
00344   //-- First can be (y,-x,0), unless direciton is (0,0,1), or close
00345   CLHEP::Hep3Vector PerpAxis1;
00346   if(fabs(fabs(_direction.z())-1.) > 0.1) {
00347     if (ALIUtils::debug >= 99) ALIUtils::dump3v( _direction, "_direction1");
00348     PerpAxis1 = CLHEP::Hep3Vector(_direction.y(),-_direction.x(),0.);
00349   } else {
00350     if (ALIUtils::debug >= 99) ALIUtils::dump3v( _direction, "_direction2");
00351     PerpAxis1 = CLHEP::Hep3Vector(_direction.z(),0.,-_direction.y());
00352   }
00353   if (ALIUtils::debug >= 4) ALIUtils::dump3v( PerpAxis1, "1st perpendicular direction of DDet");
00354 
00355   CLHEP::Hep3Vector PerpAxis2 = _direction.cross(PerpAxis1);
00356   if (ALIUtils::debug >= 4) ALIUtils::dump3v( PerpAxis2, "2nd perpendicular direction of DDet");
00357 
00358   //---------- Shift 
00359   CLHEP::Hep3Vector pointold = _point;
00360   _point += shiftAxis1*PerpAxis1;
00361   _point += shiftAxis2*PerpAxis2;
00362   if(_point != pointold && ALIUtils::debug >= 3 ) {
00363     ALIUtils::dump3v( _point-pointold, "CHANGE point");
00364   }
00365 
00366   //---------- Deviate
00367   CLHEP::Hep3Vector direcold = _direction;
00368   _direction.rotate(deviAxis1, PerpAxis1);
00369   _direction.rotate(deviAxis2, PerpAxis2);
00370   if(_direction != direcold && ALIUtils::debug >= 3) {
00371     ALIUtils::dump3v( _direction-direcold, "CHANGE direction");
00372   }
00373 
00374 } 
00375 */
00376 
00377 
00378 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00379 //@@ 
00380 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00381 void LightRay::dumpData(const ALIstring& str) const
00382 {
00383   std::cout << str << std::endl;
00384   ALIUtils::dump3v( _point, "$$ LightRay point: ");
00385   ALIUtils::dump3v( _direction, "$$ LightRay direction: ");
00386   /*
00387   CLHEP::Hep3Vector dirn = _direction;
00388   dirn.rotateZ( -23.72876*3.1415926/180.);
00389   ALIUtils::dump3v( dirn, "$$ LightRay direction: ");
00390   */
00391 }
/data/refman/pasoursint/CMSSW_4_4_5_patch3/src/Alignment/CocoaModel/src/LightRay.cc
