Line data Source code
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 22 : 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 0 : for(UInt_t i = 0; i < 3; ++i)
46 0 : v[i] = p2[i] - p1[i];
47 0 : }
48 :
49 : #if 0
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 : }
57 : #endif
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 0 : Vector(pts[0], pts[2], v[0]);
66 : //v[1] will be cross-product (v[2] x v[0]).
67 0 : Vector(pts[0], pts[1], v[2]);
68 0 : }
69 :
70 : //___________________________________________________________________
71 : Double_t Length(const Double_t (&v)[3])
72 : {
73 0 : return TMath::Sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
74 : }
75 :
76 : #if 0
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 : //______________________________________________________________________________
86 : #endif
87 : void CrossProduct(const Double_t (&v1)[3], const Double_t (&v2)[3], Double_t (&v3)[3])
88 : {
89 0 : v3[0] = v1[1] * v2[2] - v2[1] * v1[2];
90 0 : v3[1] = v1[2] * v2[0] - v2[2] * v1[0];
91 0 : v3[2] = v1[0] * v2[1] - v2[0] * v1[1];
92 0 : }
93 :
94 : //___________________________________________________________________
95 : void Normalize(Double_t (&v)[3])
96 : {
97 0 : const Double_t len = Length(v);
98 0 : if(len < 1E-10)//Threshold?
99 0 : throw std::runtime_error("Zero len vector");
100 0 : v[0] /= len;
101 0 : v[1] /= len;
102 0 : v[2] /= len;
103 0 : }
104 :
105 : //______________________________________________________________________________
106 : void Normalize(Double_t (&v)[3][3])
107 : {
108 0 : for(UInt_t i = 0; i < 3; ++i)
109 0 : Normalize(v[i]);
110 0 : }
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 0 : r[0] = v[0][0] * u[0][0] + v[1][0] * u[1][0] + v[2][0] * u[2][0];
119 0 : r[1] = v[0][0] * u[0][1] + v[1][0] * u[1][1] + v[2][0] * u[2][1];
120 0 : r[2] = v[0][0] * u[0][2] + v[1][0] * u[1][2] + v[2][0] * u[2][2];
121 :
122 0 : r[3] = v[0][1] * u[0][0] + v[1][1] * u[1][0] + v[2][1] * u[2][0];
123 0 : r[4] = v[0][1] * u[0][1] + v[1][1] * u[1][1] + v[2][1] * u[2][1];
124 0 : r[5] = v[0][1] * u[0][2] + v[1][1] * u[1][2] + v[2][1] * u[2][2];
125 :
126 0 : r[6] = v[0][2] * u[0][0] + v[1][2] * u[1][0] + v[2][2] * u[2][0];
127 0 : r[7] = v[0][2] * u[0][1] + v[1][2] * u[1][1] + v[2][2] * u[2][1];
128 0 : r[8] = v[0][2] * u[0][2] + v[1][2] * u[1][2] + v[2][2] * u[2][2];
129 0 : }
130 :
131 : //___________________________________________________________________
132 : void Rotate(const Double_t (&r)[9], const Double_t (&u)[3], Double_t (&v)[3])
133 : {
134 0 : v[0] = r[0] * u[0] + r[1] * u[1] + r[2] * u[2];
135 0 : v[1] = r[3] * u[0] + r[4] * u[1] + r[5] * u[2];
136 0 : v[2] = r[6] * u[0] + r[7] * u[1] + r[8] * u[2];
137 0 : }
138 :
139 : //___________________________________________________________________
140 : void Rotate(const Double_t (&r)[9], const Double_t (&u)[3][3], Double_t (&v)[3][3])
141 : {
142 0 : for(UInt_t i = 0; i < 3; ++i)
143 0 : Rotate(r, u[i], v[i]);
144 0 : }
145 :
146 : /*
147 : PrintVector, PrintMatrix, Translate are used in "debug" mode only.
148 : */
149 : //___________________________________________________________________
150 : void PrintVector(const Double_t (&v)[3])
151 : {
152 0 : std::cout<<v[0]<<' '<<v[1]<<' '<<v[2]<<std::endl;
153 0 : }
154 :
155 : //___________________________________________________________________
156 : void PrintMatrix(const Double_t (&u)[3][3])
157 : {
158 0 : for(UInt_t i = 0; i < 3; ++i)
159 0 : PrintVector(u[i]);
160 0 : }
161 :
162 : //___________________________________________________________________
163 : void Translate(const Double_t (&t)[3], const Double_t (&u)[3], Double_t (&v)[3])
164 : {
165 0 : for(UInt_t i = 0; i < 3; ++i)
166 0 : v[i] = u[i] + t[i];
167 0 : }
168 :
169 : //___________________________________________________________________
170 : void Translate(const Double_t (&t)[3], const Double_t (&u)[3][3], Double_t (&v)[3][3])
171 : {
172 0 : for(UInt_t i = 0; i < 3; ++i)
173 0 : Translate(t, u[i], v[i]);
174 0 : }
175 :
176 : }
177 :
178 : //______________________________________________________________________________
179 : AliPHOSModuleMisalignment::
180 : AliPHOSModuleMisalignment(const AliPHOSGeometry & geom, Bool_t debug)
181 0 : : fDebug(debug),
182 0 : fAngles(),
183 0 : fCenters(),
184 0 : fModule(),
185 0 : fU(),
186 0 : fV()
187 0 : {
188 : //Ctor.
189 : //Extract ideal module transformations from AliPHOSGeometry.
190 :
191 : //Angles.
192 0 : for (UInt_t module = 0; module < kModules; ++module)
193 0 : for (UInt_t axis = 0; axis < 3; ++axis)
194 0 : for (UInt_t angle = 0; angle < 2; ++angle)
195 0 : fAngles[module][axis][angle] = geom.GetModuleAngle(module, axis, angle);
196 : //Translations.
197 0 : for (UInt_t module = 0; module < kModules; ++module)
198 0 : for (UInt_t axis = 0; axis < 3; ++axis)
199 0 : fCenters[module][axis] = geom.GetModuleCenter(module, axis);
200 : //Points, will be rotated/translated using module angle/center.
201 0 : fModule[0][0] = -geom.GetNPhi() / 2. * geom.GetCellStep() + geom.GetCellStep() / 2.;
202 0 : fModule[0][1] = -geom.GetNZ() / 2. * geom.GetCellStep() + geom.GetCellStep() / 2.;
203 0 : 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 0 : fModule[1][0] = fModule[0][0];
206 0 : fModule[1][1] = -fModule[0][1] - geom.GetCellStep();
207 0 : fModule[1][2] = -22.61;
208 :
209 0 : fModule[2][0] = -fModule[0][0] - 7 * geom.GetCellStep();
210 0 : fModule[2][1] = fModule[0][1];
211 0 : fModule[2][2] = -22.61;
212 0 : }
213 :
214 : //______________________________________________________________________________
215 : AliPHOSModuleMisalignment::~AliPHOSModuleMisalignment()
216 0 : {
217 0 : }
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 0 : const AliSurveyPoint * pt0 = static_cast<AliSurveyPoint *>(points->FindObject(name0));
227 0 : const AliSurveyPoint * pt1 = static_cast<AliSurveyPoint *>(points->FindObject(name1));
228 0 : const AliSurveyPoint * pt2 = static_cast<AliSurveyPoint *>(points->FindObject(name2));
229 :
230 0 : if (!pt0 || !pt1 || !pt2) {
231 0 : Warning("AliPHOSModuleData::DeltaTransformation",
232 : "One of points not found in TObjArray");
233 0 : return;
234 : }
235 :
236 : //Transform fModule using angle and translation for module number "module".
237 : //fU.
238 0 : FindIdealModule(module);
239 : //Extract coordinates from survey.
240 : //fV.
241 0 : FindRealModule(pt0, pt1, pt2);
242 : //Find delta, using ideal module (fU) and survey data (fV).
243 0 : FindDelta(delta);
244 0 : }
245 :
246 : //______________________________________________________________________________
247 : void AliPHOSModuleMisalignment::FindIdealModule(UInt_t module)
248 : {
249 : //Ideal module described by fU.
250 0 : TGeoHMatrix matrix;
251 : //Rotation.
252 0 : const TGeoRotation r("",
253 0 : fAngles[module][0][0], fAngles[module][0][1],
254 0 : fAngles[module][1][0], fAngles[module][1][1],
255 0 : fAngles[module][2][0], fAngles[module][2][1]);
256 0 : matrix.SetRotation(r.GetRotationMatrix());
257 : //Translation.
258 0 : matrix.SetDx(fCenters[module][0]);
259 0 : matrix.SetDy(fCenters[module][1]);
260 0 : matrix.SetDz(fCenters[module][2]);
261 : //Find ideal module's points.
262 0 : matrix.LocalToMaster(fModule[0], fU[0]);
263 0 : matrix.LocalToMaster(fModule[1], fU[1]);
264 0 : matrix.LocalToMaster(fModule[2], fU[2]);
265 0 : }
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 0 : fV[0][0] = pt0->GetX() * scale;
278 0 : fV[0][1] = pt0->GetY() * scale;
279 0 : fV[0][2] = pt0->GetZ() * scale;
280 :
281 0 : fV[1][0] = pt1->GetX() * scale;
282 0 : fV[1][1] = pt1->GetY() * scale;
283 0 : fV[1][2] = pt1->GetZ() * scale;
284 :
285 0 : fV[2][0] = pt2->GetX() * scale;
286 0 : fV[2][1] = pt2->GetY() * scale;
287 0 : fV[2][2] = pt2->GetZ() * scale;
288 0 : }
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 0 : Double_t u[3][3] = {};
296 0 : FindVectors(fU, u);
297 : //Find cross product u2 x u0 and save it in u[2].
298 0 : 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 0 : Double_t v[3][3] = {};
304 0 : FindVectors(fV, v);
305 : //Find cross product (v2 x v0) and save it in v[2].
306 0 : 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 0 : Normalize(u);
315 0 : Normalize(v);
316 0 : } catch (const std::exception & e) {
317 : //One of lengths is zero (in principle, impossible, just to be neat).
318 0 : 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 0 : }
323 :
324 : //2. Rotation matrix.
325 0 : Double_t r[9] = {};
326 0 : FindRotation(u, v, r);
327 0 : delta->SetRotation(r);
328 :
329 : #if 1
330 :
331 : //3. Now, rotate fU and look, what translation I have to add.
332 0 : Double_t idealRotated[3] = {};
333 0 : Rotate(r, fU[0], idealRotated);
334 :
335 0 : delta->SetDx(fV[0][0] - idealRotated[0]);
336 0 : delta->SetDy(fV[0][1] - idealRotated[1]);
337 0 : delta->SetDz(fV[0][2] - idealRotated[2]);
338 :
339 0 : if (fDebug) {
340 0 : const Double_t shifts[3] =
341 0 : {fV[0][0] - idealRotated[0],
342 0 : fV[0][1] - idealRotated[1],
343 0 : fV[0][2] - idealRotated[2]};
344 :
345 0 : Double_t test1[3][3] = {};
346 0 : Rotate(r, fU, test1);
347 0 : Double_t test2[3][3] = {};
348 0 : Translate(shifts, test1, test2);
349 0 : std::cout<<"ideal:\n";
350 0 : PrintMatrix(fU);
351 0 : std::cout<<"misaligned:\n";
352 0 : PrintMatrix(test2);
353 0 : std::cout<<"real:\n";
354 0 : PrintMatrix(fV);
355 0 : }
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 0 : }
|
