CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ALILine.cc
Go to the documentation of this file.
1 // COCOA class implementation file
2 //Id: ALILine.C
3 //CAT: Fit
4 //
5 // History: v1.0
6 // Pedro Arce
7 
8 #include "Alignment/CocoaModel/interface/ALILine.h"
9 #include "Alignment/CocoaUtilities/interface/CocoaGlobals.h"
10 #include "Alignment/CocoaModel/interface/ALIPlane.h"
11 #include "Alignment/CocoaUtilities/interface/ALIUtils.h"
12 #include <cstdlib>
13 
14 
15 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
16 //@@ Constructor
17 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
18 ALILine::ALILine( const CLHEP::Hep3Vector& point, const CLHEP::Hep3Vector& direction )
19 {
20  _point = point;
21  _direction = direction* (1/direction.mag());
22  if (ALIUtils::debug >= 9) ALIUtils::dump3v( _point, " ALILine point _point = ");
23  if (ALIUtils::debug >= 9) ALIUtils::dump3v( _direction, " ALILine direction _direction = ");
24 }
25 
26 
27 
28 
29 CLHEP::Hep3Vector ALILine::intersect( const ALILine& l2, bool notParallel )
30 {
31 
32  if(ALIUtils::debug >= 3 ){
33  std::cout << "***** ALILine::intersect (constructor) two lines: " << std::endl
34  << " line1: " << *this << std::endl
35  << " line2: " << l2 << std::endl;
36  }
37  CLHEP::Hep3Vector inters;
38  //--- Check that they are not parallel
39  double acosvec = vec().dot(l2.vec()) / vec().mag() / l2.vec().mag();
40 
41  if(ALIUtils::debug >= 5 ){
42 
43  std::cout << "\t Determination of acosvec = vec().dot(l2.vec()) / vec().mag() / l2.vec().mag() "
44  << std::endl << std::endl;
45 // std::cout << "\t vec().dot = " << vec().dot << std::endl;
46  std::cout << "\t l2.vec() = " << l2.vec() << std::endl;
47 
48  std::cout << "\t vec().mag() = " << vec().mag() << std::endl;
49  std::cout << "\t l2.vec().mag() = " << l2.vec().mag() << std::endl << std::endl;
50 
51  }
52 
53  if(ALIUtils::debug >= 5 ) std::cout << " acosvec = " << acosvec << std::endl;
54  if( 1 - fabs(acosvec) < 1E-8 ) {
55  if( notParallel ) {
56  std::cerr << " !!!EXITING ALILine::intersect: two lines are parallel"
57  << std::endl;
58  exit(1);
59  } else {
60  if(ALIUtils::debug >= 5 ) std::cout << " !!! ALILine::intersect: two lines are parallel (no errors)" << std::endl;
61  //gcc2952 inters = CLHEP::Hep3Vector( DBL_MAX, DBL_MAX, DBL_MAX );
62  inters = CLHEP::Hep3Vector( ALI_DBL_MAX, ALI_DBL_MAX, ALI_DBL_MAX );
63  }
64  } else {
65 
66 
67  // **************************************************************** //
68  // **************************************************************** //
69  // **************************************************************** //
70  // Determination of Fact //
71  // //
72  // Problem : 3D quantity was determined by doing calculation with //
73  // the 2D projections of the std::vectors. It is possible //
74  // for projection in a particular plane to be 0 //
75  // //
76  // Solution : Test for problem and redo calculation if necessary by //
77  // projecting into a different plane //
78  // **************************************************************** //
79  // **************************************************************** //
80  // **************************************************************** //
81 
82 
83 
84  ALIdouble fact = ( vec().y() * l2.pt().x() - vec().x() * l2.pt().y()
85  - vec().y() * pt().x() + vec().x() * pt().y() )
86  / ( vec().x() * l2.vec().y() - vec().y() * l2.vec().x() );
87 
88 
89 
90  if(ALIUtils::debug >= 2 ){
91  std::cout << std::endl << std::endl << "*** START CALC OF FACT ***" << std::endl;
92  std::cout << " ==================" << std::endl << std::endl;
93  std::cout << "*** Determination of fact ->";
94  std::cout << "\t fact = (" << vec().y() * l2.pt().x() - vec().x() * l2.pt().y()
95  - vec().y() * pt().x() + vec().x() * pt().y() << "/";
96 
97  std::cout << vec().x() * l2.vec().y() - vec().y() * l2.vec().x() << ") = " << fact << std::endl;
98 
99  }
100 
101  // ALIdouble old_fact = fact;
102  ALIdouble old_fact_denominator = vec().x() * l2.vec().y() - vec().y() * l2.vec().x();
103  //- ALIdouble old_fact2 = 0.;
104  ALIdouble old_fact_denominator2 = 999;
105 
106 
107  if(fabs(fact) > 1e8 || fabs(old_fact_denominator) < 1.e-10)
108  {
109 
110  // do calculation by rotating 90 degrees into xy plane
111  // Z-> X
112  // X-> -Z
113 
114  if(ALIUtils::debug >= 2 && fabs(old_fact_denominator) < 1.e-10) std::cout << " ** Problem: old_fact_denominator -> " << old_fact_denominator << std::endl;
115 
116  if(ALIUtils::debug >= 2 && fabs(fact) > 1e8) std::cout << " ** Problem: fact -> " << fact << std::endl;
117 
118  if(ALIUtils::debug >= 2 && (vec().x() * l2.vec().y() - vec().y() * l2.vec().x()) == 0) std::cout << " ** Division by 0 !!! " << std::endl;
119  if(ALIUtils::debug >= 2 ) std::cout << " ** Must rotate to yz plane for calculation (X-> Z) ";
120 
121 
122  fact = ( -1 * vec().y() * l2.pt().z() + vec().z() * l2.pt().y()
123  + vec().y() * pt().z() - vec().z() * pt().y() )
124  / ( -1*vec().z() * l2.vec().y() + vec().y() * l2.vec().z() );
125 
126  if(ALIUtils::debug >= 2 ) std::cout << "\t -- 1st Recalculation of fact in yz plane = " << fact << std::endl;
127 
128  old_fact_denominator2 = -1*vec().z() * l2.vec().y() + vec().y() * l2.vec().z();
129 
130 
131  if(fabs(-1*vec().z() * l2.vec().y() + vec().y() * l2.vec().z())< 1.e-10)
132  {
133 
134  if(ALIUtils::debug >= 2 ) std::cout << " ** Must rotate to xz plane for calculation (Y-> -Z) ";
135  if(ALIUtils::debug >= 2 && fabs(old_fact_denominator2) < 1.e-10) std::cout << " ** Problem: old_fact_denominator2 -> " << old_fact_denominator2 << std::endl;
136  //- if(ALIUtils::debug >= 2 && fabs(old_fact2) > 1.e8) std::cout << " ** Problem: old_fact2 -> " << old_fact2 << std::endl;
137 
138  fact = ( -1*vec().z() * l2.pt().x() + vec().x() * l2.pt().z()
139  + vec().z() * pt().x() - vec().x() * pt().z() )
140  / ( -1*vec().x() * l2.vec().z() + vec().z() * l2.vec().x() );
141 
142  if(ALIUtils::debug >= 2 ) std::cout << "\t -- 2nd Recalculation of fact in xz plane = " << fact << std::endl;
143 
144 
145  }
146  else {
147  if(ALIUtils::debug >= 2 ) std::cout << "*!* 2nd calculation sufficient" << std::endl;
148  }
149 
150  }
151  else
152  {
153 
154  if(ALIUtils::debug >= 2 ) std::cout << "*!* Standard calculation - things are fine" << std::endl;
155  if(ALIUtils::debug >= 2 ) std::cout << "\t ----> fact = ( " << vec().y() * l2.pt().x() - vec().x() *
156  l2.pt().y() - vec().y() * pt().x() + vec().x() * pt().y() << " / " << vec().x() * l2.vec().y()
157  - vec().y()* l2.vec().x() << " ) = " << fact << std::endl;
158  }
159 
160 
161 
162  // **************************************************************** //
163  // **************************************************************** //
164  // **************************************************************** //
165  // Finished With Fact //
166  // **************************************************************** //
167  // **************************************************************** //
168  // **************************************************************** //
169 
170 
171  inters = l2.pt() + fact * l2.vec();
172 
173  if(ALIUtils::debug >= 2 ){
174  std::cout << "Determination of intersection = l2.pt() + fact * l2.vec()" << std::endl << std::endl;
175  ALIUtils::dump3v( l2.pt(), "\t --> l2.pt() = ");
176  std::cout << "\t --> fact = " << fact << std::endl;
177  std::cout << "\t --> l2.vec() = " << l2.vec() << std::endl;
178  ALIUtils::dump3v( inters, "***\t --> ALILine::intersection at: ");
179  }
180 
181  }
182 
183  return inters;
184 
185 }
186 
187 std::ostream& operator<<(std::ostream& out, const ALILine& li)
188 {
189  out << " ALILine point " << li._point << std::endl;
190  out << " ALILine direc " << li._direction;
191 
192  return out;
193 }
194 
195 
196 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
197 //@@ intersect: Intersect a LightRay with a plane and change thePoint to the intersection point
198 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
199 
200 CLHEP::Hep3Vector ALILine::intersect( const ALIPlane& plane, bool notParallel)
201 {
202  if(ALIUtils::debug >= 5) std::cout << "***** ALILine::intersect WITH PLANE" << std::endl;
203  if(ALIUtils::debug >= 4) {
204  ALIUtils::dump3v( plane.point(), "plane point");
205  ALIUtils::dump3v( plane.normal(), "plane normal");
206  }
207 
208  //---------- Check that they intersect
209  if( fabs( plane.normal()*_direction ) < 1.E-10 ) {
210  std::cerr << " !!!! INTERSECTION NOT POSSIBLE: LightRay is perpendicular to plane " << std::endl;
211  std::cerr << " plane.normal()*direction() = " << plane.normal()*_direction << std::endl;
212  ALIUtils::dump3v( _direction, "LightRay direction ");
213  ALIUtils::dump3v( plane.normal(), "plane normal ");
214  exit(1);
215  }
216 
217  //---------- Get intersection point between LightRay and plane
218  //- if(ALIUtils::debug >= 4) std::cout << " ALILine::intersect WITH LightRay" << std::endl;
219 
220  CLHEP::Hep3Vector vtemp = plane.point() - _point;
221  if(ALIUtils::debug >= 4) ALIUtils::dump3v( vtemp, "n_r = point - point_plane");
222 
223  ALIdouble dtemp = _direction * plane.normal();
224  if(ALIUtils::debug >= 4) std::cout << " lightray* plate normal" << dtemp << std::endl;
225 
226  if ( dtemp != 0. ) {
227  dtemp = (vtemp * plane.normal()) / dtemp;
228  if(ALIUtils::debug >= 4) std::cout << " n_r*plate normal (dtemp) : " << dtemp << std::endl;
229  if(ALIUtils::debug >= 4) std::cout << " Old vtemp : " << vtemp << std::endl;
230  } else std::cerr << "!!! LightRay: Intersect With Plane: plane and light ray parallel: no intersection" << std::endl;
231 
232  vtemp = _direction * dtemp;
233  if(ALIUtils::debug >= 4) ALIUtils::dump3v( vtemp, "n_r scaled (vtemp) : ");
234  if(ALIUtils::debug >= 4) ALIUtils::dump3v( CLHEP::Hep3Vector(dtemp), "dtemp analog to vtemp : ");
235 
236  CLHEP::Hep3Vector inters = vtemp + _point;
237  if(ALIUtils::debug >= 4) ALIUtils::dump3v( inters, "intersection point = vtemp + _point");
238  if(ALIUtils::debug >= 4) ALIUtils::dump3v( vtemp, "vtemp = ");
239  if(ALIUtils::debug >= 4) ALIUtils::dump3v( _point, "_point = ");
240 
241  return inters;
242 }
243 
ALIdouble
long double ALIdouble
Definition: CocoaGlobals.h:11
ALILine::_point
CLHEP::Hep3Vector _point
Definition: ALILine.h:33
ALIUtils::debug
static ALIint debug
Definition: ALIUtils.h:35
operator<<
std::ostream & operator<<(std::ostream &out, const ALILine &li)
Definition: ALILine.cc:187
detailsBasic3DVector::z
float float float z
Definition: extBasic3DVector.h:15
ecal_dqm_sourceclient-live_cfg.cerr
tuple cerr
Definition: ecal_dqm_sourceclient-live_cfg.py:89
ALILine::ALILine
ALILine()
Definition: ALILine.h:19
ALILine::pt
const CLHEP::Hep3Vector & pt() const
Definition: ALILine.h:27
ALIPlane::point
const CLHEP::Hep3Vector & point() const
Definition: ALIPlane.h:23
dbtoconf.out
tuple out
Definition: dbtoconf.py:99
fact
const double fact
Definition: NuclearInteractionFTFSimulator.cc:61
ALILine::_direction
CLHEP::Hep3Vector _direction
Definition: ALILine.h:34
ALIUtils::dump3v
static void dump3v(const CLHEP::Hep3Vector &vec, const std::string &msg)
Definition: ALIUtils.cc:61
ALILine::intersect
CLHEP::Hep3Vector intersect(const ALILine &l2, bool notParallel=0)
Definition: ALILine.cc:29
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
x
Definition: DDAxes.h:10
ALIPlane::normal
const CLHEP::Hep3Vector & normal() const
Definition: ALIPlane.h:24
point
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
ALILine::vec
const CLHEP::Hep3Vector & vec() const
Definition: ALILine.h:28
ALI_DBL_MAX
const double ALI_DBL_MAX
Definition: CocoaGlobals.h:24
Generated for CMSSW Reference Manual by doxygen 1.8.5