fixes
[u/mrichter/AliRoot.git] / FMD / AliFMDSurveyToAlignObjs.cxx
CommitLineData
09b6c804 1//
2// Class to take survey data and
3// transform that to alignment objects.
4//
5// FMD
6//
b2e6f0b0 7#include "AliFMDSurveyToAlignObjs.h"
f567c3ce 8#include "AliLog.h"
b2e6f0b0 9#include "AliSurveyPoint.h"
10#include <TGraph2DErrors.h>
11#include <TF2.h>
12#include <TVector3.h>
13#include <iostream>
14#include <iomanip>
15#include <TMath.h>
16#include <TRotation.h>
17#include <TGeoMatrix.h>
18#include <TGeoManager.h>
19#include <TGeoPhysicalNode.h>
20#include "AliFMDGeometry.h"
21
22//____________________________________________________________________
faf80567 23Double_t
24AliFMDSurveyToAlignObjs::GetUnitFactor() const
25{
26 // Returns the conversion factor from the measured values to
27 // centimeters.
f567c3ce 28 if (!fSurveyObj) return 0;
faf80567 29 TString units(fSurveyObj->GetUnits());
30 if (units.CompareTo("mm", TString::kIgnoreCase) == 0) return .1;
31 else if (units.CompareTo("cm", TString::kIgnoreCase) == 0) return 1.;
32 else if (units.CompareTo("m", TString::kIgnoreCase) == 0) return 100.;
33 return 1;
34}
35
36//____________________________________________________________________
b2e6f0b0 37Bool_t
faf80567 38AliFMDSurveyToAlignObjs::GetPoint(const char* name,
39 TVector3& point,
40 TVector3& error) const
b2e6f0b0 41{
faf80567 42 // Get named point. On return, point will contain the point
43 // coordinates in centimeters, and error will contain the
44 // meassurement errors in centimeters too. If no point is found,
45 // returns false, otherwise true.
f567c3ce 46 if (!fSurveyPoints) return kFALSE;
47
faf80567 48 Double_t unit = GetUnitFactor();
f567c3ce 49 if (unit == 0) return kFALSE;
50
faf80567 51 TObject* obj = fSurveyPoints->FindObject(name);
52 if (!obj) return kFALSE;
53
f567c3ce 54 AliSurveyPoint* p = static_cast<AliSurveyPoint*>(obj);
55 point.SetXYZ(unit * p->GetX(),
56 unit * p->GetY(),
57 unit * p->GetZ());
58 error.SetXYZ(unit * p->GetPrecisionX(),
59 unit * p->GetPrecisionY(),
60 unit * p->GetPrecisionZ());
61 return kTRUE;
faf80567 62}
b2e6f0b0 63
faf80567 64//____________________________________________________________________
65Bool_t
66AliFMDSurveyToAlignObjs::CalculatePlane(const TVector3& a,
67 const TVector3& b,
68 const TVector3& c,
f567c3ce 69 Double_t depth,
faf80567 70 Double_t* trans,
71 Double_t* rot) const
72{
09b6c804 73 //
74 // Calculate the plane translation and rotation from 3 survey points
75 //
76 // Parameters:
77 // a 1st Survey point
78 // b 2nd Survey point
79 // c 3rd Survey point
80 // trans Translation vector
81 // rot Rotation matrix (direction cosines)
82 //
83 // Return:
84 //
85 //
86
faf80567 87 // Vector a->b, b->c, and normal to plane defined by these two
88 // vectors.
89 TVector3 ab(b-a), bc(c-b);
b2e6f0b0 90
faf80567 91 // Normal vector to the plane of the fiducial marks obtained
92 // as cross product of the two vectors on the plane d0^d1
93 TVector3 nn(ab.Cross(bc));
94 if (nn.Mag() < 1e-8) {
f567c3ce 95 Info("CalculatePlane", "Normal vector is null vector");
faf80567 96 return kFALSE;
97 }
b2e6f0b0 98
faf80567 99 // We express the plane in Hessian normal form.
100 //
101 // n x = -p,
102 //
103 // where n is the normalised normal vector given by
104 //
105 // n_x = a / l, n_y = b / l, n_z = c / l, p = d / l
106 //
107 // with l = sqrt(a^2+b^2+c^2) and a, b, c, and d are from the
108 // normal plane equation
109 //
110 // ax + by + cz + d = 0
111 //
112 // Normalize
113 TVector3 n(nn.Unit());
f567c3ce 114 // Double_t p = - (n * a);
faf80567 115
116 // The center of the square with the fiducial marks as the
117 // corners. The mid-point of one diagonal - md. Used to get the
118 // center of the surveyd box.
119 TVector3 md(a + c);
120 md *= 1/2.;
f567c3ce 121
faf80567 122 // The center of the box.
f567c3ce 123 TVector3 orig(md - depth * n);
faf80567 124 trans[0] = orig[0];
125 trans[1] = orig[1];
126 trans[2] = orig[2];
127
128 // Normalize the spanning vectors
129 TVector3 uab(ab.Unit());
130 TVector3 ubc(bc.Unit());
131
132 for (size_t i = 0; i < 3; i++) {
f567c3ce 133 rot[i * 3 + 0] = ubc[i];
134 rot[i * 3 + 1] = uab[i];
faf80567 135 rot[i * 3 + 2] = n[i];
136 }
137 return kTRUE;
138}
139
faf80567 140//____________________________________________________________________
f567c3ce 141Bool_t
142AliFMDSurveyToAlignObjs::FitPlane(const TObjArray& points,
143 const TObjArray& errors,
144 Double_t /* depth */,
145 Double_t* trans,
146 Double_t* rot) const
faf80567 147{
09b6c804 148 //
149 // Calculate the plane rotation and translation by doing a fit of
150 // the plane equation to the surveyed points. At least 4 points
151 // must be passed in the @a points array with corresponding errors
152 // in the array @a errors. The arrays are assumed to contain
153 // TVector3 objects.
154 //
155 // Parameters:
156 // points Array surveyed positions
157 // errors Array of errors corresponding to @a points
158 // depth Survey targets depth (perpendicular to the plane)
159 // trans On return, translation of the plane
160 // rot On return, rotation (direction cosines) of the plane
161 //
162 // Return:
163 // @c true on success, @c false otherwise
164 //
165
f567c3ce 166 Int_t nPoints = points.GetEntries();
167 if (nPoints < 4) {
168 AliError(Form("Cannot fit a plane equation to less than 4 survey points, "
169 "got only %d", nPoints));
faf80567 170 return kFALSE;
171 }
f567c3ce 172
faf80567 173 TGraph2DErrors g;
b2e6f0b0 174 // Loop and fill graph
f567c3ce 175 for (int i = 0; i < nPoints; i++) {
176 TVector3* p = static_cast<TVector3*>(points.At(i));
177 TVector3* e = static_cast<TVector3*>(errors.At(i));
b2e6f0b0 178
f567c3ce 179 if (!p || !e) continue;
180
181 g.SetPoint(i, p->X(), p->Y(), p->Z());
182 g.SetPointError(i, e->X(), e->Y(), e->Z());
183 }
184
b2e6f0b0 185 // Check that we have enough points
186 if (g.GetN() < 4) {
f567c3ce 187 AliError(Form("Only got %d survey points - no good for plane fit",
188 g.GetN()));
b2e6f0b0 189 return kFALSE;
190 }
191
192 // Next, declare fitting function and fit to graph.
193 // Fit to the plane equation:
194 //
195 // ax + by + cz + d = 0
196 //
197 // or
198 //
199 // z = - ax/c - by/c - d/c
200 //
f567c3ce 201 TF2 f("plane", "-[0]*x-[1]*y-[2]",
202 g.GetXmin(), g.GetXmax(), g.GetYmin(), g.GetYmax());
b2e6f0b0 203 g.Fit(&f, "Q");
f567c3ce 204
b2e6f0b0 205 // Now, extract the normal and offset
206 TVector3 nv(f.GetParameter(0), f.GetParameter(1), 1);
207 TVector3 n(nv.Unit());
208 Double_t p = -f.GetParameter(2);
f567c3ce 209
b2e6f0b0 210 // Create two vectors spanning the plane
211 TVector3 a(1, 0, f.Eval(1, 0)-p);
212 TVector3 b(0, -1, f.Eval(0, -1)-p);
213 TVector3 ua(a.Unit());
214 TVector3 ub(b.Unit());
f567c3ce 215 // Double_t angAb = ua.Angle(ub);
faf80567 216 // PrintVector("ua: ", ua);
217 // PrintVector("ub: ", ub);
218 // std::cout << "Angle: " << angAb * 180 / TMath::Pi() << std::endl;
f567c3ce 219
b2e6f0b0 220 for (size_t i = 0; i < 3; i++) {
221 rot[i * 3 + 0] = ua[i];
222 rot[i * 3 + 1] = ub[i];
223 rot[i * 3 + 2] = n[i];
224 }
f567c3ce 225
b2e6f0b0 226 // The intersection of the plane is given by (0, 0, -d/c)
227 trans[0] = 0;
228 trans[1] = 0;
229 trans[2] = p;
f567c3ce 230
b2e6f0b0 231 return kTRUE;
232}
233
f567c3ce 234
b2e6f0b0 235//____________________________________________________________________
236Bool_t
f567c3ce 237AliFMDSurveyToAlignObjs::MakeDelta(const char* path,
09b6c804 238 const Double_t* rot,
239 const Double_t* trans,
f567c3ce 240 TGeoHMatrix& delta) const
b2e6f0b0 241{
09b6c804 242 //
243 // Create a delta transform from a global rotation matrix and
244 // translation.
245 //
246 // Parameters:
247 // path Path of element to transform.
248 // rot Rotation matrix (direction cosines)
249 // trans Translation
250 // delta On return, the delta transform
251 //
252 // Return:
253 // Newly
254 //
f567c3ce 255 if (!gGeoManager) return kFALSE;
256 if (!gGeoManager->cd(path)) return kFALSE;
b2e6f0b0 257
b2e6f0b0 258
f567c3ce 259 TGeoMatrix* global = gGeoManager->GetCurrentMatrix();
260#if 0
261 PrintRotation(Form("%s rot:", global->GetName()),global->GetRotationMatrix());
262 PrintVector(Form("%s trans:", global->GetName()),global->GetTranslation());
263#endif
b2e6f0b0 264
f567c3ce 265 return MakeDelta(global, rot, trans, delta);
266}
267
268//____________________________________________________________________
269Bool_t
09b6c804 270AliFMDSurveyToAlignObjs::MakeDelta(const TGeoMatrix* global,
271 const Double_t* rot,
272 const Double_t* trans,
f567c3ce 273 TGeoHMatrix& delta) const
274{
09b6c804 275 //
276 // Create a delta transform from a global rotation matrix and
277 // translation.
278 //
279 // Parameters:
280 // global Global matrix of element to transform.
281 // rot Rotation matrix (direction cosines)
282 // trans Translation
283 // delta On return, the delta transform
284 //
285 // Return:
286 // Newly
287 //
b2e6f0b0 288 TGeoHMatrix* geoM = new TGeoHMatrix;
b2e6f0b0 289 geoM->SetTranslation(trans);
f567c3ce 290 geoM->SetRotation(rot);
291
292 delta = global->Inverse();
293 delta.MultiplyLeft(geoM);
294
295 return true;
296}
297
298//____________________________________________________________________
299Bool_t
300AliFMDSurveyToAlignObjs::GetFMD1Plane(Double_t* rot, Double_t* trans) const
301{
09b6c804 302 //
303 // Get the FMD1 plane from the survey points
304 //
305 // Parameters:
306 // rot Rotation matrix (direction cosines)
307 // trans Translation
308 //
309 // Return:
310 // @c true on success, @c false otherwise.
311 //
f567c3ce 312
313 // The possile survey points
314 TVector3 icb, ict, ocb, oct, dummy;
315 Int_t missing = 0;
316 if (!GetPoint("V0L_ICB", icb, dummy)) missing++;
317 if (!GetPoint("V0L_ICT", ict, dummy)) missing++;
318 if (!GetPoint("V0L_OCB", ocb, dummy)) missing++;
319 if (!GetPoint("V0L_OCT", oct, dummy)) missing++;
320
321 // Check that we have enough points
322 if (missing > 1) {
323 AliWarning(Form("Only got %d survey points - no good for FMD1 plane",
324 4-missing));
325 return kFALSE;
326 }
327
328#if 0
329 const char* lidN = "FMD1_lid_mat0";
330 TGeoMatrix* lidM = static_cast<TGeoMatrix*>(gGeoManager->GetListOfMatrices()
331 ->FindObject(lidN));
332 if (!lidM) {
333 AliError(Form("Couldn't find FMD1 lid transformation %s", lidN));
334 return kFALSE;
335 }
336
337 const Double_t* lidT = lidM->GetTranslation();
338 Double_t lidZ = lidT[2];
339 Double_t off = lidZ-3.3;
340#else
341 Double_t off = 0;
b2e6f0b0 342#endif
f567c3ce 343
344 if (!CalculatePlane(ocb, icb, ict, off, trans, rot)) return kFALSE;
345 // PrintRotation("FMD1 rotation:", rot);
346 // PrintVector("FMD1 translation:", trans);
347
348 return kTRUE;
349}
350
351//____________________________________________________________________
352Bool_t
353AliFMDSurveyToAlignObjs::DoFMD1()
354{
09b6c804 355 //
356 // Do the FMD1 analysis. We have 4 survey targets on V0-A on the
357 // C-side. These are
358 //
359 // - V0A_ICT In-side, C-side, top.
360 // - V0A_ICB In-side, C-side, bottom.
361 // - V0A_OCT Out-side, C-side, top.
362 // - V0A_OCB Out-side, C-side, bottom.
363 //
364 // These 4 survey targets sit 3.3mm over the V0-A C-side surface, or
365 // 3.3mm over the back surface of FMD1.
366 //
367 // Since these are really sitting on a plane, we can use the method
368 // proposed by the CORE offline.
369 //
370 // Return:
371 // @c true on success, @c false otherwise.
372 //
373
f567c3ce 374 // Do the FMD1 stuff
375 Double_t rot[9], trans[3];
376 if (!GetFMD1Plane(rot, trans)) return kFALSE;
377 // const char* path = "/ALIC_1/F1MT_1/FMD1_lid_0";
378
379 // TGeoHMatrix delta;
380 TGeoTranslation global(0,0,324.670);
381 if (!MakeDelta(&global, rot, trans, fFMD1Delta))
382 return kFALSE;
b2e6f0b0 383
f567c3ce 384 // PrintRotation("FMD1 delta rotation:", fFMD1Delta.GetRotationMatrix());
385 // PrintVector("FMD1 delta translation:", fFMD1Delta.GetTranslation());
386
387 return kTRUE;
388}
389
390//____________________________________________________________________
391Bool_t
392AliFMDSurveyToAlignObjs::GetFMD2Plane(Double_t* rot, Double_t* trans) const
393{
09b6c804 394 //
395 // Get the surveyed plane corresponding to the backside of FMD2.
396 // The plane is done as a best fit of the plane equation to at least
397 // 4 of the available survey points.
398 //
399 // Parameters:
400 // rot Rotation matrix (direction cosines)
401 // trans Translation vector.
402 //
403 // Return:
404 // @c true on success, @c false otherwise
405 //
f567c3ce 406
407 // The possible survey points
408 const char* names[] = { "FMD2_ITOP", "FMD2_OTOP",
409 "FMD2_IBOTM", "FMD2_OBOTM",
410 "FMD2_IBOT", "FMD2_OBOT",
411 0 };
412 const char** name = names;
413
414 TObjArray points;
415 TObjArray errors;
b2e6f0b0 416
f567c3ce 417 // Loop and fill graph
418 while (*name) {
419 TVector3 p, e;
420 if (!GetPoint(*name++, p, e)) continue;
421
422 points.Add(new TVector3(p));
423 errors.Add(new TVector3(e));
424 }
425 if (points.GetEntries() < 4) {
426 AliWarning(Form("Only got %d survey points - no good for FMD2 plane",
427 points.GetEntries()));
428 return kFALSE;
429 }
430
431 return FitPlane(points, errors, 0, trans, rot);
432}
433
434#define M(I,J) rot[(J-1) * 3 + (I-1)]
435//____________________________________________________________________
436Bool_t
437AliFMDSurveyToAlignObjs::DoFMD2()
438{
09b6c804 439 //
440 // Do the FMD2 calculations. We have 6 survey points of which only
441 // 5 are normally surveyed. These are all sittings
442 //
443 // - FMD2_ITOP - In-side, top
444 // - FMD2_IBOTM - In-side, middle bottom
445 // - FMD2_IBOT - In-side, bottom
446 // - FMD2_OTOP - Out-side, top
447 // - FMD2_OBOTM - Out-side, middle bottom
448 // - FMD2_OBOT - Out-side, bottom
449 //
450 // The nominal coordinates of these retro-fitted survey stickers
451 // isn't known. Also, these stickers are put on a thin (0.3mm
452 // thick) carbon cover which flexes quite easily. This means, that
453 // to rotations and xy-translation obtained from the survey data
454 // cannot be used, and left is only the z-translation.
455 //
456 // Further more, since FMD2 to is attached to the ITS SPD thermal
457 // screen, it is questionable if the FMD2 survey will ever be used.
458 //
459 // Return:
460 // @c true on success, @c false otherwise.
461 //
462
f567c3ce 463 // Do the FMD2 stuff
464 Double_t rot[9], trans[3];
465 if (!GetFMD2Plane(rot, trans)) return kFALSE;
466 // PrintRotation("FMD2 rotation:", rot);
467 // PrintVector("FMD2 translation:", trans);
468
469 for (int i = 0; i < 3; i++) {
470 for (int j = 0; j < 3; j++) {
471 rot[i*3+j] = (i == j ? 1 : 0);
472 }
473 }
474 trans[0] = trans[1] = 0;
475 trans[2] += 0.015;
476 // PrintRotation("FMD2 rotation:", rot);
477 // PrintVector("FMD2 translation:", trans);
b2e6f0b0 478
f567c3ce 479 // TGeoHMatrix delta;
480 if (!MakeDelta("/ALIC_1/F2MT_2/FMD2_support_0/FMD2_back_cover_2",
481 rot, trans, fFMD2Delta)) return kFALSE;
b2e6f0b0 482
f567c3ce 483 // PrintRotation("FMD2 delta rotation:", fFMD2Delta.GetRotationMatrix());
484 // PrintVector("FMD2 delta translation:", fFMD2Delta.GetTranslation());
485
b2e6f0b0 486 return kTRUE;
487}
488
489//____________________________________________________________________
490void
491AliFMDSurveyToAlignObjs::Run()
492{
09b6c804 493 //
494 // Run the task.
495 //
496 //
497
b2e6f0b0 498 AliFMDGeometry* geom = AliFMDGeometry::Instance();
499 geom->Init();
500 geom->InitTransformations();
501
f567c3ce 502 DoFMD1();
b2e6f0b0 503 DoFMD2();
504}
505
506//____________________________________________________________________
f567c3ce 507Bool_t
508AliFMDSurveyToAlignObjs::CreateAlignObjs()
509{
09b6c804 510 //
511 //
512 // Method to create the alignment objects
513 //
514 // Return:
515 // @c true on success, @c false otherwise
516 //
f567c3ce 517 TClonesArray& array = *fAlignObjArray;
518 Int_t n = array.GetEntriesFast();
519
520 if (!fFMD1Delta.IsIdentity()) {
521 new (array[n++]) AliAlignObjParams("FMD1/FMD1_T", 0, fFMD1Delta, kTRUE);
522 new (array[n++]) AliAlignObjParams("FMD1/FMD1_B", 0, fFMD1Delta, kTRUE);
523 }
524 if (!fFMD2Delta.IsIdentity()) {
525 new (array[n++]) AliAlignObjParams("FMD2/FMD2_T", 0, fFMD2Delta, kTRUE);
526 new (array[n++]) AliAlignObjParams("FMD2/FMD2_B", 0, fFMD2Delta, kTRUE);
527 }
528 // array.Print();
529
530 return kTRUE;
531}
532
533//____________________________________________________________________
b2e6f0b0 534void
535AliFMDSurveyToAlignObjs::PrintVector(const char* text, const TVector3& v)
536{
09b6c804 537 //
538 // Service member function to print a vector
539 //
540 // Parameters:
541 // text Prefix text
542 // v Vector
543 //
b2e6f0b0 544 Double_t va[] = { v.X(), v.Y(), v.Z() };
545 PrintVector(text, va);
546}
547//____________________________________________________________________
548void
549AliFMDSurveyToAlignObjs::PrintVector(const char* text, const Double_t* v)
550{
09b6c804 551 //
552 // Service member function to print a vector
553 //
554 // Parameters:
555 // text Prefix text
556 // v Vector (array of 3 doubles)
557 //
b2e6f0b0 558 std::cout << text
559 << std::setw(15) << v[0]
560 << std::setw(15) << v[1]
561 << std::setw(15) << v[2]
562 << std::endl;
563}
564
565
566//____________________________________________________________________
567void
568AliFMDSurveyToAlignObjs::PrintRotation(const char* text, const Double_t* rot)
569{
09b6c804 570 //
571 // Service member function to print a rotation matrix
572 //
573 // Parameters:
574 // text Prefix text
575 // v Matrix (array of 9 doubles)
576 //
577
b2e6f0b0 578 std::cout << text << std::endl;
579 for (size_t i = 0; i < 3; i++) {
580 for (size_t j = 0; j < 3; j++)
581 std::cout << std::setw(15) << rot[i * 3 + j];
582 std::cout << std::endl;
583 }
584}
585
586//____________________________________________________________________
587//
588// EOF
589//