]>
Commit | Line | Data |
---|---|---|
d7c519b2 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* $Id: $ */ | |
17 | ||
18 | #include <stdexcept> | |
19 | #include <iostream> | |
20 | ||
21 | #include <TGeoMatrix.h> | |
22 | #include <TObjArray.h> | |
23 | #include <TString.h> | |
24 | #include <TError.h> | |
25 | #include <TMath.h> | |
26 | ||
27 | #include "AliSurveyPoint.h" | |
28 | ||
29 | #include "AliPHOSModuleMisalignment.h" | |
30 | #include "AliPHOSGeometry.h" | |
31 | ||
32 | ClassImp(AliPHOSModuleMisalignment) | |
33 | ||
34 | namespace { | |
35 | ||
36 | /* | |
37 | Tiny vector/matrix utility stuff. Operates on doubles directly. | |
38 | Instead of points and vectors I use arrays of doubles with size 3. | |
39 | To make this explicit - I use references to arrays. | |
40 | */ | |
41 | ||
42 | //___________________________________________________________________ | |
43 | void Vector(const Double_t (&p1)[3], const Double_t (&p2)[3], Double_t (&v)[3]) | |
44 | { | |
45 | for(UInt_t i = 0; i < 3; ++i) | |
46 | v[i] = p2[i] - p1[i]; | |
47 | } | |
48 | ||
76664146 | 49 | #if 0 |
d7c519b2 | 50 | //___________________________________________________________________ |
51 | void MultVector(Double_t (&v)[3], Double_t m) | |
52 | { | |
53 | v[0] *= m; | |
54 | v[1] *= m; | |
55 | v[2] *= m; | |
56 | } | |
76664146 | 57 | #endif |
d7c519b2 | 58 | |
59 | /* | |
60 | Using points name0, name1, name2 find two orthogonal vectors. | |
61 | */ | |
62 | //___________________________________________________________________ | |
63 | void FindVectors(const Double_t (&pts)[3][3], Double_t (&v)[3][3]) | |
64 | { | |
65 | Vector(pts[0], pts[2], v[0]); | |
66 | //v[1] will be cross-product (v[2] x v[0]). | |
67 | Vector(pts[0], pts[1], v[2]); | |
68 | } | |
69 | ||
70 | //___________________________________________________________________ | |
71 | Double_t Length(const Double_t (&v)[3]) | |
72 | { | |
73 | return TMath::Sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); | |
74 | } | |
75 | ||
76664146 | 76 | #if 0 |
d7c519b2 | 77 | //___________________________________________________________________ |
78 | Double_t Distance(const Double_t (&p1)[3], const Double_t (&p2)[3]) | |
79 | { | |
80 | return TMath::Sqrt((p2[0] - p1[0]) * (p2[0] - p1[0]) + | |
81 | (p2[1] - p1[1]) * (p2[1] - p1[1]) + | |
82 | (p2[2] - p1[2]) * (p2[2] - p1[2])); | |
83 | } | |
84 | ||
85 | //______________________________________________________________________________ | |
76664146 | 86 | #endif |
d7c519b2 | 87 | void CrossProduct(const Double_t (&v1)[3], const Double_t (&v2)[3], Double_t (&v3)[3]) |
88 | { | |
89 | v3[0] = v1[1] * v2[2] - v2[1] * v1[2]; | |
90 | v3[1] = v1[2] * v2[0] - v2[2] * v1[0]; | |
91 | v3[2] = v1[0] * v2[1] - v2[0] * v1[1]; | |
92 | } | |
93 | ||
94 | //___________________________________________________________________ | |
95 | void Normalize(Double_t (&v)[3]) | |
96 | { | |
97 | const Double_t len = Length(v); | |
98 | if(len < 1E-10)//Threshold? | |
99 | throw std::runtime_error("Zero len vector"); | |
100 | v[0] /= len; | |
101 | v[1] /= len; | |
102 | v[2] /= len; | |
103 | } | |
104 | ||
105 | //______________________________________________________________________________ | |
106 | void Normalize(Double_t (&v)[3][3]) | |
107 | { | |
108 | for(UInt_t i = 0; i < 3; ++i) | |
109 | Normalize(v[i]); | |
110 | } | |
111 | ||
112 | ||
113 | //___________________________________________________________________ | |
114 | void FindRotation(const Double_t (&u)[3][3], const Double_t (&v)[3][3], Double_t (&r)[9]) | |
115 | { | |
116 | //I have orthogonal vectors and very nice rotation matrix. | |
117 | //V = R * U, R = V * U ^ t | |
118 | r[0] = v[0][0] * u[0][0] + v[1][0] * u[1][0] + v[2][0] * u[2][0]; | |
119 | r[1] = v[0][0] * u[0][1] + v[1][0] * u[1][1] + v[2][0] * u[2][1]; | |
120 | r[2] = v[0][0] * u[0][2] + v[1][0] * u[1][2] + v[2][0] * u[2][2]; | |
121 | ||
122 | r[3] = v[0][1] * u[0][0] + v[1][1] * u[1][0] + v[2][1] * u[2][0]; | |
123 | r[4] = v[0][1] * u[0][1] + v[1][1] * u[1][1] + v[2][1] * u[2][1]; | |
124 | r[5] = v[0][1] * u[0][2] + v[1][1] * u[1][2] + v[2][1] * u[2][2]; | |
125 | ||
126 | r[6] = v[0][2] * u[0][0] + v[1][2] * u[1][0] + v[2][2] * u[2][0]; | |
127 | r[7] = v[0][2] * u[0][1] + v[1][2] * u[1][1] + v[2][2] * u[2][1]; | |
128 | r[8] = v[0][2] * u[0][2] + v[1][2] * u[1][2] + v[2][2] * u[2][2]; | |
129 | } | |
130 | ||
131 | //___________________________________________________________________ | |
132 | void Rotate(const Double_t (&r)[9], const Double_t (&u)[3], Double_t (&v)[3]) | |
133 | { | |
134 | v[0] = r[0] * u[0] + r[1] * u[1] + r[2] * u[2]; | |
135 | v[1] = r[3] * u[0] + r[4] * u[1] + r[5] * u[2]; | |
136 | v[2] = r[6] * u[0] + r[7] * u[1] + r[8] * u[2]; | |
137 | } | |
138 | ||
139 | //___________________________________________________________________ | |
140 | void Rotate(const Double_t (&r)[9], const Double_t (&u)[3][3], Double_t (&v)[3][3]) | |
141 | { | |
142 | for(UInt_t i = 0; i < 3; ++i) | |
143 | Rotate(r, u[i], v[i]); | |
144 | } | |
145 | ||
146 | /* | |
147 | PrintVector, PrintMatrix, Translate are used in "debug" mode only. | |
148 | */ | |
149 | //___________________________________________________________________ | |
150 | void PrintVector(const Double_t (&v)[3]) | |
151 | { | |
152 | std::cout<<v[0]<<' '<<v[1]<<' '<<v[2]<<std::endl; | |
153 | } | |
154 | ||
155 | //___________________________________________________________________ | |
156 | void PrintMatrix(const Double_t (&u)[3][3]) | |
157 | { | |
158 | for(UInt_t i = 0; i < 3; ++i) | |
159 | PrintVector(u[i]); | |
160 | } | |
161 | ||
162 | //___________________________________________________________________ | |
163 | void Translate(const Double_t (&t)[3], const Double_t (&u)[3], Double_t (&v)[3]) | |
164 | { | |
165 | for(UInt_t i = 0; i < 3; ++i) | |
166 | v[i] = u[i] + t[i]; | |
167 | } | |
168 | ||
169 | //___________________________________________________________________ | |
170 | void Translate(const Double_t (&t)[3], const Double_t (&u)[3][3], Double_t (&v)[3][3]) | |
171 | { | |
172 | for(UInt_t i = 0; i < 3; ++i) | |
173 | Translate(t, u[i], v[i]); | |
174 | } | |
175 | ||
176 | } | |
177 | ||
178 | //______________________________________________________________________________ | |
179 | AliPHOSModuleMisalignment:: | |
180 | AliPHOSModuleMisalignment(const AliPHOSGeometry & geom, Bool_t debug) | |
181 | : fDebug(debug), | |
182 | fAngles(), | |
183 | fCenters(), | |
184 | fModule(), | |
185 | fU(), | |
186 | fV() | |
187 | { | |
188 | //Ctor. | |
189 | //Extract ideal module transformations from AliPHOSGeometry. | |
190 | ||
191 | //Angles. | |
192 | for (UInt_t module = 0; module < kModules; ++module) | |
193 | for (UInt_t axis = 0; axis < 3; ++axis) | |
194 | for (UInt_t angle = 0; angle < 2; ++angle) | |
195 | fAngles[module][axis][angle] = geom.GetModuleAngle(module, axis, angle); | |
196 | //Translations. | |
197 | for (UInt_t module = 0; module < kModules; ++module) | |
198 | for (UInt_t axis = 0; axis < 3; ++axis) | |
199 | fCenters[module][axis] = geom.GetModuleCenter(module, axis); | |
200 | //Points, will be rotated/translated using module angle/center. | |
201 | fModule[0][0] = -geom.GetNPhi() / 2. * geom.GetCellStep() + geom.GetCellStep() / 2.; | |
202 | fModule[0][1] = -geom.GetNZ() / 2. * geom.GetCellStep() + geom.GetCellStep() / 2.; | |
203 | fModule[0][2] = -22.61;//This number is hardcoded, AliPHOSGeometry does not have it, | |
204 | //only 460. but this is result of transformation applied already. | |
205 | fModule[1][0] = fModule[0][0]; | |
206 | fModule[1][1] = -fModule[0][1] - geom.GetCellStep(); | |
207 | fModule[1][2] = -22.61; | |
208 | ||
209 | fModule[2][0] = -fModule[0][0] - 7 * geom.GetCellStep(); | |
210 | fModule[2][1] = fModule[0][1]; | |
211 | fModule[2][2] = -22.61; | |
212 | } | |
213 | ||
214 | //______________________________________________________________________________ | |
215 | AliPHOSModuleMisalignment::~AliPHOSModuleMisalignment() | |
216 | { | |
217 | } | |
218 | ||
219 | //______________________________________________________________________________ | |
220 | void AliPHOSModuleMisalignment:: | |
221 | DeltaTransformation(UInt_t module, const TObjArray * points, | |
222 | const TString & name0, const TString & name1, | |
223 | const TString & name2, TGeoHMatrix * delta) | |
224 | { | |
225 | //Find delta transformation to misalign module. Global transformation. | |
226 | const AliSurveyPoint * pt0 = static_cast<AliSurveyPoint *>(points->FindObject(name0)); | |
227 | const AliSurveyPoint * pt1 = static_cast<AliSurveyPoint *>(points->FindObject(name1)); | |
228 | const AliSurveyPoint * pt2 = static_cast<AliSurveyPoint *>(points->FindObject(name2)); | |
229 | ||
230 | if (!pt0 || !pt1 || !pt2) { | |
231 | Warning("AliPHOSModuleData::DeltaTransformation", | |
232 | "One of points not found in TObjArray"); | |
233 | return; | |
234 | } | |
235 | ||
236 | //Transform fModule using angle and translation for module number "module". | |
237 | //fU. | |
238 | FindIdealModule(module); | |
239 | //Extract coordinates from survey. | |
240 | //fV. | |
241 | FindRealModule(pt0, pt1, pt2); | |
242 | //Find delta, using ideal module (fU) and survey data (fV). | |
243 | FindDelta(delta); | |
244 | } | |
245 | ||
246 | //______________________________________________________________________________ | |
247 | void AliPHOSModuleMisalignment::FindIdealModule(UInt_t module) | |
248 | { | |
249 | //Ideal module described by fU. | |
250 | TGeoHMatrix matrix; | |
251 | //Rotation. | |
252 | const TGeoRotation r("", | |
253 | fAngles[module][0][0], fAngles[module][0][1], | |
254 | fAngles[module][1][0], fAngles[module][1][1], | |
255 | fAngles[module][2][0], fAngles[module][2][1]); | |
256 | matrix.SetRotation(r.GetRotationMatrix()); | |
257 | //Translation. | |
258 | matrix.SetDx(fCenters[module][0]); | |
259 | matrix.SetDy(fCenters[module][1]); | |
260 | matrix.SetDz(fCenters[module][2]); | |
261 | //Find ideal module's points. | |
262 | matrix.LocalToMaster(fModule[0], fU[0]); | |
263 | matrix.LocalToMaster(fModule[1], fU[1]); | |
264 | matrix.LocalToMaster(fModule[2], fU[2]); | |
265 | } | |
266 | ||
267 | //______________________________________________________________________________ | |
268 | void AliPHOSModuleMisalignment::FindRealModule(const AliSurveyPoint * pt0, | |
269 | const AliSurveyPoint * pt1, | |
270 | const AliSurveyPoint * pt2) | |
271 | { | |
272 | //Real module - fV. | |
273 | //Survey is in millimeters. | |
274 | //AliPHOSGeometry is in centimeters. | |
275 | const Double_t scale = 0.1; | |
276 | ||
277 | fV[0][0] = pt0->GetX() * scale; | |
278 | fV[0][1] = pt0->GetY() * scale; | |
279 | fV[0][2] = pt0->GetZ() * scale; | |
280 | ||
281 | fV[1][0] = pt1->GetX() * scale; | |
282 | fV[1][1] = pt1->GetY() * scale; | |
283 | fV[1][2] = pt1->GetZ() * scale; | |
284 | ||
285 | fV[2][0] = pt2->GetX() * scale; | |
286 | fV[2][1] = pt2->GetY() * scale; | |
287 | fV[2][2] = pt2->GetZ() * scale; | |
288 | } | |
289 | ||
290 | //______________________________________________________________________________ | |
291 | void AliPHOSModuleMisalignment::FindDelta(TGeoHMatrix * delta)const | |
292 | { | |
293 | //Find rotation and translation wich can | |
294 | //convert fU into fV (ideal module points into points from survey). | |
295 | Double_t u[3][3] = {}; | |
296 | FindVectors(fU, u); | |
297 | //Find cross product u2 x u0 and save it in u[2]. | |
298 | CrossProduct(u[2], u[0], u[1]); | |
299 | /* | |
300 | const Double_t lenXideal = Length(u[0]); | |
301 | const Double_t lenZideal = Length(u[2]); | |
302 | */ | |
303 | Double_t v[3][3] = {}; | |
304 | FindVectors(fV, v); | |
305 | //Find cross product (v2 x v0) and save it in v[2]. | |
306 | CrossProduct(v[2], v[0], v[1]); | |
307 | /* | |
308 | const Double_t lenXreal = Length(v[0]); | |
309 | const Double_t lenZreal = Length(v[2]); | |
310 | */ | |
311 | //Now, find rotation matrix. | |
312 | //1. Normalize vectors in u and v. | |
313 | try { | |
314 | Normalize(u); | |
315 | Normalize(v); | |
316 | } catch (const std::exception & e) { | |
317 | //One of lengths is zero (in principle, impossible, just to be neat). | |
318 | Error("AliPHOSModuleMisalignment::FindDelta", | |
319 | "\tone of vectors from ideal or real\n\tpoints have zero size\n" | |
320 | "\tzero misalignment will be created"); | |
321 | return; | |
322 | } | |
323 | ||
324 | //2. Rotation matrix. | |
325 | Double_t r[9] = {}; | |
326 | FindRotation(u, v, r); | |
327 | delta->SetRotation(r); | |
328 | ||
329 | #if 1 | |
330 | ||
331 | //3. Now, rotate fU and look, what translation I have to add. | |
332 | Double_t idealRotated[3] = {}; | |
333 | Rotate(r, fU[0], idealRotated); | |
334 | ||
335 | delta->SetDx(fV[0][0] - idealRotated[0]); | |
336 | delta->SetDy(fV[0][1] - idealRotated[1]); | |
337 | delta->SetDz(fV[0][2] - idealRotated[2]); | |
338 | ||
339 | if (fDebug) { | |
340 | const Double_t shifts[3] = | |
341 | {fV[0][0] - idealRotated[0], | |
342 | fV[0][1] - idealRotated[1], | |
343 | fV[0][2] - idealRotated[2]}; | |
344 | ||
345 | Double_t test1[3][3] = {}; | |
346 | Rotate(r, fU, test1); | |
347 | Double_t test2[3][3] = {}; | |
348 | Translate(shifts, test1, test2); | |
349 | std::cout<<"ideal:\n"; | |
350 | PrintMatrix(fU); | |
351 | std::cout<<"misaligned:\n"; | |
352 | PrintMatrix(test2); | |
353 | std::cout<<"real:\n"; | |
354 | PrintMatrix(fV); | |
355 | } | |
356 | ||
357 | #endif | |
358 | ||
359 | #if 0 | |
360 | //3. Now, rotate fU and look, what translation I have to add. | |
361 | Double_t idealRotated[3][3] = {}; | |
362 | Rotate(r, fU, idealRotated); | |
363 | //Because of measurement errors, distances | |
364 | //between points has errors. I can try to split | |
365 | //this difference (and make "final errors" smaller). | |
366 | Double_t zShift[3] = {}; | |
367 | Vector(fV[0], fV[1], zShift); | |
368 | Normalize(zShift); | |
369 | ||
370 | Double_t xShift[3] = {}; | |
371 | Vector(fV[0], fV[2], xShift); | |
372 | Normalize(xShift); | |
373 | ||
374 | MultVector(zShift, 0.5 * (lenZreal - lenZideal)); | |
375 | MultVector(xShift, 0.5 * (lenXreal - lenXideal)); | |
376 | ||
377 | Double_t pt1[3] = {}; | |
378 | Translate(zShift, fV[0], pt1); | |
379 | Double_t pt2[3] = {}; | |
380 | Translate(xShift, pt1, pt2); | |
381 | ||
382 | Double_t shifts[] = {pt2[0] - idealRotated[0][0], | |
383 | pt2[1] - idealRotated[0][1], | |
384 | pt2[2] - idealRotated[0][2]}; | |
385 | ||
386 | delta->SetDx(shifts[0]); | |
387 | delta->SetDy(shifts[1]); | |
388 | delta->SetDz(shifts[2]); | |
389 | ||
390 | if (fDebug) { | |
391 | Double_t idealTr[3][3] = {}; | |
392 | Translate(shifts, idealRotated, idealTr); | |
393 | ||
394 | std::cout<<"misaligned:\n"; | |
395 | PrintMatrix(idealTr); | |
396 | std::cout<<"ideal1 "<<Distance(idealTr[0], idealTr[1])<<std::endl; | |
397 | std::cout<<"ideal2 "<<Distance(idealTr[0], idealTr[2])<<std::endl; | |
398 | std::cout<<"real:\n"; | |
399 | PrintMatrix(fV); | |
400 | std::cout<<"real1 "<<Distance(fV[0], fV[1])<<std::endl; | |
401 | std::cout<<"real2 "<<Distance(fV[0], fV[2])<<std::endl; | |
402 | } | |
403 | #endif | |
404 | } |