2 // Class to take survey data and
3 // transform that to alignment objects.
7 #include "AliFMDSurveyToAlignObjs.h"
9 #include "AliSurveyPoint.h"
10 #include <TGraph2DErrors.h>
16 #include <TRotation.h>
17 #include <TGeoMatrix.h>
18 #include <TGeoManager.h>
19 #include <TGeoPhysicalNode.h>
20 #include "AliFMDGeometry.h"
22 //____________________________________________________________________
24 AliFMDSurveyToAlignObjs::GetUnitFactor() const
26 // Returns the conversion factor from the measured values to
28 if (!fSurveyObj) return 0;
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.;
36 //____________________________________________________________________
38 AliFMDSurveyToAlignObjs::GetPoint(const char* name,
40 TVector3& error) const
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.
46 if (!fSurveyPoints) return kFALSE;
48 Double_t unit = GetUnitFactor();
49 if (unit == 0) return kFALSE;
51 TObject* obj = fSurveyPoints->FindObject(name);
52 if (!obj) return kFALSE;
54 AliSurveyPoint* p = static_cast<AliSurveyPoint*>(obj);
55 point.SetXYZ(unit * p->GetX(),
58 error.SetXYZ(unit * p->GetPrecisionX(),
59 unit * p->GetPrecisionY(),
60 unit * p->GetPrecisionZ());
64 //____________________________________________________________________
66 AliFMDSurveyToAlignObjs::CalculatePlane(const TVector3& a,
74 // Calculate the plane translation and rotation from 3 survey points
80 // trans Translation vector
81 // rot Rotation matrix (direction cosines)
87 // Vector a->b, b->c, and normal to plane defined by these two
89 TVector3 ab(b-a), bc(c-b);
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) {
95 Info("CalculatePlane", "Normal vector is null vector");
99 // We express the plane in Hessian normal form.
103 // where n is the normalised normal vector given by
105 // n_x = a / l, n_y = b / l, n_z = c / l, p = d / l
107 // with l = sqrt(a^2+b^2+c^2) and a, b, c, and d are from the
108 // normal plane equation
110 // ax + by + cz + d = 0
113 TVector3 n(nn.Unit());
114 // Double_t p = - (n * a);
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.
122 // The center of the box.
123 TVector3 orig(md - depth * n);
128 // Normalize the spanning vectors
129 TVector3 uab(ab.Unit());
130 TVector3 ubc(bc.Unit());
132 for (size_t i = 0; i < 3; i++) {
133 rot[i * 3 + 0] = ubc[i];
134 rot[i * 3 + 1] = uab[i];
135 rot[i * 3 + 2] = n[i];
140 //____________________________________________________________________
142 AliFMDSurveyToAlignObjs::FitPlane(const TObjArray& points,
143 const TObjArray& errors,
144 Double_t /* depth */,
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
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
163 // @c true on success, @c false otherwise
166 Int_t nPoints = points.GetEntries();
168 AliError(Form("Cannot fit a plane equation to less than 4 survey points, "
169 "got only %d", nPoints));
174 // Loop and fill graph
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));
179 if (!p || !e) continue;
181 g.SetPoint(i, p->X(), p->Y(), p->Z());
182 g.SetPointError(i, e->X(), e->Y(), e->Z());
185 // Check that we have enough points
187 AliError(Form("Only got %d survey points - no good for plane fit",
192 // Next, declare fitting function and fit to graph.
193 // Fit to the plane equation:
195 // ax + by + cz + d = 0
199 // z = - ax/c - by/c - d/c
201 TF2 f("plane", "-[0]*x-[1]*y-[2]",
202 g.GetXmin(), g.GetXmax(), g.GetYmin(), g.GetYmax());
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);
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());
215 // Double_t angAb = ua.Angle(ub);
216 // PrintVector("ua: ", ua);
217 // PrintVector("ub: ", ub);
218 // std::cout << "Angle: " << angAb * 180 / TMath::Pi() << std::endl;
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];
226 // The intersection of the plane is given by (0, 0, -d/c)
235 //____________________________________________________________________
237 AliFMDSurveyToAlignObjs::MakeDelta(const char* path,
239 const Double_t* trans,
240 TGeoHMatrix& delta) const
243 // Create a delta transform from a global rotation matrix and
247 // path Path of element to transform.
248 // rot Rotation matrix (direction cosines)
250 // delta On return, the delta transform
255 if (!gGeoManager) return kFALSE;
256 if (!gGeoManager->cd(path)) return kFALSE;
259 TGeoMatrix* global = gGeoManager->GetCurrentMatrix();
261 PrintRotation(Form("%s rot:", global->GetName()),global->GetRotationMatrix());
262 PrintVector(Form("%s trans:", global->GetName()),global->GetTranslation());
265 return MakeDelta(global, rot, trans, delta);
268 //____________________________________________________________________
270 AliFMDSurveyToAlignObjs::MakeDelta(const TGeoMatrix* global,
272 const Double_t* trans,
273 TGeoHMatrix& delta) const
276 // Create a delta transform from a global rotation matrix and
280 // global Global matrix of element to transform.
281 // rot Rotation matrix (direction cosines)
283 // delta On return, the delta transform
288 TGeoHMatrix* geoM = new TGeoHMatrix;
289 geoM->SetTranslation(trans);
290 geoM->SetRotation(rot);
292 delta = global->Inverse();
293 delta.MultiplyLeft(geoM);
298 //____________________________________________________________________
300 AliFMDSurveyToAlignObjs::GetFMD1Plane(Double_t* rot, Double_t* trans) const
303 // Get the FMD1 plane from the survey points
306 // rot Rotation matrix (direction cosines)
310 // @c true on success, @c false otherwise.
313 // The possile survey points
314 TVector3 icb, ict, ocb, oct, dummy;
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++;
321 // Check that we have enough points
323 AliWarning(Form("Only got %d survey points - no good for FMD1 plane",
329 const char* lidN = "FMD1_lid_mat0";
330 TGeoMatrix* lidM = static_cast<TGeoMatrix*>(gGeoManager->GetListOfMatrices()
333 AliError(Form("Couldn't find FMD1 lid transformation %s", lidN));
337 const Double_t* lidT = lidM->GetTranslation();
338 Double_t lidZ = lidT[2];
339 Double_t off = lidZ-3.3;
344 if (!CalculatePlane(ocb, icb, ict, off, trans, rot)) return kFALSE;
345 // PrintRotation("FMD1 rotation:", rot);
346 // PrintVector("FMD1 translation:", trans);
351 //____________________________________________________________________
353 AliFMDSurveyToAlignObjs::DoFMD1()
356 // Do the FMD1 analysis. We have 4 survey targets on V0-A on the
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.
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.
367 // Since these are really sitting on a plane, we can use the method
368 // proposed by the CORE offline.
371 // @c true on success, @c false otherwise.
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";
379 // TGeoHMatrix delta;
380 TGeoTranslation global(0,0,324.670);
381 if (!MakeDelta(&global, rot, trans, fFMD1Delta))
384 // PrintRotation("FMD1 delta rotation:", fFMD1Delta.GetRotationMatrix());
385 // PrintVector("FMD1 delta translation:", fFMD1Delta.GetTranslation());
390 //____________________________________________________________________
392 AliFMDSurveyToAlignObjs::GetFMD2Plane(Double_t* rot, Double_t* trans) const
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.
400 // rot Rotation matrix (direction cosines)
401 // trans Translation vector.
404 // @c true on success, @c false otherwise
407 // The possible survey points
408 const char* names[] = { "FMD2_ITOP", "FMD2_OTOP",
409 "FMD2_IBOTM", "FMD2_OBOTM",
410 "FMD2_IBOT", "FMD2_OBOT",
412 const char** name = names;
417 // Loop and fill graph
420 if (!GetPoint(*name++, p, e)) continue;
422 points.Add(new TVector3(p));
423 errors.Add(new TVector3(e));
425 if (points.GetEntries() < 4) {
426 AliWarning(Form("Only got %d survey points - no good for FMD2 plane",
427 points.GetEntries()));
431 return FitPlane(points, errors, 0, trans, rot);
434 #define M(I,J) rot[(J-1) * 3 + (I-1)]
435 //____________________________________________________________________
437 AliFMDSurveyToAlignObjs::DoFMD2()
440 // Do the FMD2 calculations. We have 6 survey points of which only
441 // 5 are normally surveyed. These are all sittings
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
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.
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.
460 // @c true on success, @c false otherwise.
464 Double_t rot[9], trans[3];
465 if (!GetFMD2Plane(rot, trans)) return kFALSE;
466 // PrintRotation("FMD2 rotation:", rot);
467 // PrintVector("FMD2 translation:", trans);
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);
474 trans[0] = trans[1] = 0;
476 // PrintRotation("FMD2 rotation:", rot);
477 // PrintVector("FMD2 translation:", trans);
479 // TGeoHMatrix delta;
480 if (!MakeDelta("/ALIC_1/F2MT_2/FMD2_support_0/FMD2_back_cover_2",
481 rot, trans, fFMD2Delta)) return kFALSE;
483 // PrintRotation("FMD2 delta rotation:", fFMD2Delta.GetRotationMatrix());
484 // PrintVector("FMD2 delta translation:", fFMD2Delta.GetTranslation());
489 //____________________________________________________________________
491 AliFMDSurveyToAlignObjs::Run()
498 AliFMDGeometry* geom = AliFMDGeometry::Instance();
500 geom->InitTransformations();
506 //____________________________________________________________________
508 AliFMDSurveyToAlignObjs::CreateAlignObjs()
512 // Method to create the alignment objects
515 // @c true on success, @c false otherwise
517 TClonesArray& array = *fAlignObjArray;
518 Int_t n = array.GetEntriesFast();
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);
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);
533 //____________________________________________________________________
535 AliFMDSurveyToAlignObjs::PrintVector(const char* text, const TVector3& v)
538 // Service member function to print a vector
544 Double_t va[] = { v.X(), v.Y(), v.Z() };
545 PrintVector(text, va);
547 //____________________________________________________________________
549 AliFMDSurveyToAlignObjs::PrintVector(const char* text, const Double_t* v)
552 // Service member function to print a vector
556 // v Vector (array of 3 doubles)
559 << std::setw(15) << v[0]
560 << std::setw(15) << v[1]
561 << std::setw(15) << v[2]
566 //____________________________________________________________________
568 AliFMDSurveyToAlignObjs::PrintRotation(const char* text, const Double_t* rot)
571 // Service member function to print a rotation matrix
575 // v Matrix (array of 9 doubles)
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;
586 //____________________________________________________________________