Reference Manual

Go to the documentation of this file.
00001 //   COCOA class implementation file
00002 //Id:  ALILine.C
00003 //CAT: Fit
00004 //
00005 //   History: v1.0 
00006 //   Pedro Arce
00007 
00008 #include "Alignment/CocoaModel/interface/ALILine.h"
00009 #include "Alignment/CocoaUtilities/interface/CocoaGlobals.h" 
00010 #include "Alignment/CocoaModel/interface/ALIPlane.h"
00011 #include "Alignment/CocoaUtilities/interface/ALIUtils.h"
00012 #include <cstdlib>
00013 
00014 
00015 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00016 //@@ Constructor
00017 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00018 ALILine::ALILine( const CLHEP::Hep3Vector& point, const CLHEP::Hep3Vector& direction ) 
00019 {
00020   _point = point;
00021   _direction = direction* (1/direction.mag());
00022   if (ALIUtils::debug >= 9) ALIUtils::dump3v( _point, " ALILine point _point = ");
00023   if (ALIUtils::debug >= 9) ALIUtils::dump3v( _direction, " ALILine direction _direction = ");
00024 }
00025 
00026 
00027 
00028 
00029 CLHEP::Hep3Vector ALILine::intersect( const ALILine& l2, bool notParallel ) 
00030 { 
00031   
00032   if(ALIUtils::debug >= 3 ){
00033     std::cout << "***** ALILine::intersect (constructor) two lines: " << std::endl 
00034                         << " line1: " << *this << std::endl
00035                         << " line2: " << l2 << std::endl;
00036   }
00037   CLHEP::Hep3Vector inters; 
00038   //--- Check that they are not parallel 
00039   double acosvec = vec().dot(l2.vec()) / vec().mag() / l2.vec().mag();
00040   
00041   if(ALIUtils::debug >= 5 ){
00042   
00043         std::cout << "\t Determination of acosvec = vec().dot(l2.vec()) / vec().mag() / l2.vec().mag() "
00044         << std::endl << std::endl;
00045 //      std::cout << "\t vec().dot = " << vec().dot << std::endl; 
00046         std::cout << "\t l2.vec() = " << l2.vec() << std::endl;
00047 
00048         std::cout << "\t vec().mag() = " << vec().mag() << std::endl;
00049         std::cout << "\t l2.vec().mag() = " << l2.vec().mag() << std::endl << std::endl;
00050   
00051                         }
00052   
00053   if(ALIUtils::debug >= 5 ) std::cout << " acosvec = " << acosvec << std::endl;
00054   if( 1 - fabs(acosvec) < 1E-8 ) {
00055     if( notParallel ) {
00056       std::cerr << " !!!EXITING ALILine::intersect: two lines are parallel"
00057         << std::endl;
00058       exit(1);
00059     } else {
00060       if(ALIUtils::debug >= 5 ) std::cout << " !!! ALILine::intersect: two lines are parallel (no errors)" << std::endl;
00061       //gcc2952      inters = CLHEP::Hep3Vector( DBL_MAX, DBL_MAX, DBL_MAX );
00062       inters = CLHEP::Hep3Vector( ALI_DBL_MAX, ALI_DBL_MAX, ALI_DBL_MAX );
00063     }
00064   } else { 
00065 
00066     
00067  //     ****************************************************************        //  
00068  //     ****************************************************************        //
00069  //     ****************************************************************        //
00070  //                          Determination of Fact                              //
00071  //                                                                             //
00072  //     Problem :  3D quantity was determined by doing calculation with         //
00073  //                the 2D projections of the std::vectors.  It is possible              //
00074  //                for projection in a particular plane to be 0                 //
00075  //                                                                             //
00076  //     Solution : Test for problem and redo calculation if necessary by        //
00077  //                projecting into a different plane                            //
00078  //     ****************************************************************        //  
00079  //     ****************************************************************        //
00080  //     ****************************************************************        //
00081   
00082  
00083         
00084     ALIdouble fact = ( vec().y() * l2.pt().x() - vec().x() * l2.pt().y()
00085        - vec().y() * pt().x() + vec().x() * pt().y() )
00086       / ( vec().x() * l2.vec().y() - vec().y() * l2.vec().x() );
00087 
00088 
00089 
00090   if(ALIUtils::debug >= 2 ){
00091        std::cout << std::endl << std::endl << "*** START CALC OF FACT ***" << std::endl;
00092        std::cout << "    ==================" << std::endl << std::endl;            
00093        std::cout << "*** Determination of fact ->";      
00094        std::cout << "\t fact = (" << vec().y() * l2.pt().x() - vec().x() * l2.pt().y()
00095        - vec().y() * pt().x() + vec().x() * pt().y() << "/";
00096        
00097        std::cout <<  vec().x() * l2.vec().y() - vec().y() * l2.vec().x() << ") = " << fact << std::endl;
00098  
00099                                 }
00100         
00101   //    ALIdouble old_fact = fact;
00102     ALIdouble old_fact_denominator = vec().x() * l2.vec().y() - vec().y() * l2.vec().x();                               
00103     //-    ALIdouble old_fact2 = 0.;
00104     ALIdouble old_fact_denominator2 = 999;
00105                
00106 
00107    if(fabs(fact) > 1e8 || fabs(old_fact_denominator) < 1.e-10)
00108     {
00109       
00110       // do calculation by rotating 90 degrees into xy plane
00111       // Z-> X
00112       // X-> -Z
00113       
00114        if(ALIUtils::debug >= 2 && fabs(old_fact_denominator) < 1.e-10) std::cout << " ** Problem:  old_fact_denominator -> " << old_fact_denominator << std::endl;
00115        
00116        if(ALIUtils::debug >= 2 && fabs(fact) > 1e8) std::cout << " ** Problem: fact -> " << fact << std::endl;
00117        
00118        if(ALIUtils::debug >= 2 && (vec().x() * l2.vec().y() - vec().y() * l2.vec().x()) == 0) std::cout << " ** Division by 0 !!! " << std::endl;       
00119        if(ALIUtils::debug >= 2 ) std::cout << " ** Must rotate to yz plane for calculation (X-> Z) ";
00120  
00121     
00122         fact = ( -1 * vec().y() * l2.pt().z() + vec().z() * l2.pt().y()
00123         + vec().y() * pt().z() - vec().z() * pt().y() )
00124         / ( -1*vec().z() * l2.vec().y() + vec().y() * l2.vec().z() );
00125       
00126         if(ALIUtils::debug >= 2 ) std::cout << "\t -- 1st Recalculation of fact in yz plane = " << fact << std::endl;  
00127 
00128         old_fact_denominator2 = -1*vec().z() * l2.vec().y() + vec().y() * l2.vec().z();
00129 
00130         
00131         if(fabs(-1*vec().z() * l2.vec().y() + vec().y() * l2.vec().z())< 1.e-10)
00132         {
00133           
00134           if(ALIUtils::debug >= 2 ) std::cout << " ** Must rotate to xz plane for calculation (Y-> -Z) ";
00135             if(ALIUtils::debug >= 2 && fabs(old_fact_denominator2) < 1.e-10) std::cout << " ** Problem: old_fact_denominator2 -> " << old_fact_denominator2 << std::endl;       
00136           //-              if(ALIUtils::debug >= 2 && fabs(old_fact2) > 1.e8) std::cout << " ** Problem: old_fact2 -> " << old_fact2 << std::endl;
00137           
00138           fact = ( -1*vec().z() * l2.pt().x() + vec().x() * l2.pt().z()
00139                    + vec().z() * pt().x() - vec().x() * pt().z() )
00140             / ( -1*vec().x() * l2.vec().z() + vec().z() * l2.vec().x() );
00141           
00142            if(ALIUtils::debug >= 2 ) std::cout << "\t -- 2nd Recalculation of fact in xz plane = " << fact << std::endl;  
00143 
00144 
00145         }
00146         else {
00147           if(ALIUtils::debug >= 2 ) std::cout << "*!* 2nd calculation sufficient" << std::endl;
00148         }
00149        
00150     }    
00151     else
00152     {
00153 
00154      if(ALIUtils::debug >= 2 ) std::cout << "*!* Standard calculation - things are fine" << std::endl;
00155      if(ALIUtils::debug >= 2 ) std::cout << "\t ----> fact = ( " << vec().y() * l2.pt().x() - vec().x() *
00156      l2.pt().y() - vec().y() * pt().x() + vec().x() * pt().y() << " / " << vec().x() * l2.vec().y() 
00157      - vec().y()* l2.vec().x() << " ) = " << fact << std::endl;    
00158     }
00159 
00160           
00161           
00162  //     ****************************************************************        //  
00163  //     ****************************************************************        //
00164  //     ****************************************************************        //
00165  //                             Finished With Fact                                      //
00166  //     ****************************************************************        //  
00167  //     ****************************************************************        //
00168  //     ****************************************************************        //        
00169           
00170           
00171      inters =  l2.pt() + fact * l2.vec();
00172      
00173      if(ALIUtils::debug >= 2 ){           
00174         std::cout << "Determination of intersection = l2.pt() + fact * l2.vec()" << std::endl << std::endl;
00175         ALIUtils::dump3v( l2.pt(), "\t --> l2.pt() = ");
00176         std::cout << "\t --> fact = " << fact << std::endl;
00177         std::cout << "\t --> l2.vec() = " << l2.vec() << std::endl;
00178         ALIUtils::dump3v( inters, "***\t --> ALILine::intersection at: ");
00179      }
00180 
00181   }
00182  
00183   return inters;
00184  
00185 }
00186 
00187 std::ostream& operator<<(std::ostream& out, const ALILine& li)
00188 {
00189   out << " ALILine point " << li._point << std::endl;
00190   out << " ALILine direc " << li._direction;
00191 
00192   return out;  
00193 }
00194 
00195 
00196 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00197 //@@ intersect: Intersect a LightRay with a plane and change thePoint to the intersection point
00198 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00199 
00200 CLHEP::Hep3Vector ALILine::intersect( const ALIPlane& plane, bool notParallel)
00201 {
00202   if(ALIUtils::debug >= 5) std::cout << "***** ALILine::intersect WITH PLANE" << std::endl;
00203   if(ALIUtils::debug >= 4) {
00204     ALIUtils::dump3v( plane.point(), "plane point"); 
00205     ALIUtils::dump3v( plane.normal(), "plane normal"); 
00206   }
00207   
00208   //---------- Check that they intersect
00209   if( fabs( plane.normal()*_direction ) < 1.E-10 ) {
00210     std::cerr << " !!!! INTERSECTION NOT POSSIBLE: LightRay is perpendicular to plane " << std::endl;
00211     std::cerr << " plane.normal()*direction() = " << plane.normal()*_direction << std::endl;
00212     ALIUtils::dump3v( _direction, "LightRay direction ");
00213     ALIUtils::dump3v( plane.normal(), "plane normal ");
00214     exit(1);
00215   }
00216 
00217   //---------- Get intersection point between LightRay and plane
00218   //-   if(ALIUtils::debug >= 4) std::cout << " ALILine::intersect WITH LightRay" << std::endl;
00219         
00220   CLHEP::Hep3Vector vtemp = plane.point() - _point;
00221         if(ALIUtils::debug >= 4) ALIUtils::dump3v( vtemp, "n_r = point  - point_plane"); 
00222         
00223   ALIdouble dtemp = _direction * plane.normal();
00224         if(ALIUtils::debug >= 4) std::cout << " lightray* plate normal" << dtemp << std::endl;
00225         
00226   if ( dtemp != 0. ) {
00227       dtemp = (vtemp * plane.normal()) / dtemp;
00228       if(ALIUtils::debug >= 4)  std::cout << " n_r*plate normal (dtemp) : " << dtemp << std::endl;
00229       if(ALIUtils::debug >= 4)  std::cout << " Old vtemp : " << vtemp << std::endl;
00230         } else std::cerr << "!!! LightRay: Intersect With Plane: plane and light ray parallel: no intersection" << std::endl;     
00231                 
00232   vtemp = _direction * dtemp;
00233         if(ALIUtils::debug >= 4) ALIUtils::dump3v( vtemp, "n_r scaled (vtemp) : "); 
00234         if(ALIUtils::debug >= 4) ALIUtils::dump3v( dtemp, "dtemp analog to vtemp : "); 
00235     
00236   CLHEP::Hep3Vector inters = vtemp + _point;
00237         if(ALIUtils::debug >= 4) ALIUtils::dump3v( inters, "intersection point = vtemp + _point"); 
00238         if(ALIUtils::debug >= 4) ALIUtils::dump3v( vtemp, "vtemp =  "); 
00239         if(ALIUtils::debug >= 4) ALIUtils::dump3v( _point, "_point =  "); 
00240 
00241   return inters;
00242 }
00243
/data/refman/pasoursint/CMSSW_5_3_0/src/Alignment/CocoaModel/src/ALILine.cc
