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