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