Removing useless flag.
[u/mrichter/AliRoot.git] / PHOS / AliPHOSModuleMisalignment.cxx
CommitLineData
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
32ClassImp(AliPHOSModuleMisalignment)
33
34namespace {
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//______________________________________________________________________________
179AliPHOSModuleMisalignment::
180AliPHOSModuleMisalignment(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//______________________________________________________________________________
215AliPHOSModuleMisalignment::~AliPHOSModuleMisalignment()
216{
217}
218
219//______________________________________________________________________________
220void AliPHOSModuleMisalignment::
221DeltaTransformation(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//______________________________________________________________________________
247void 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//______________________________________________________________________________
268void 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//______________________________________________________________________________
291void 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}