1 #include "AliFMDSurveyToAlignObjs.h"
3 #include "AliSurveyPoint.h"
4 #include <TGraph2DErrors.h>
10 #include <TRotation.h>
11 #include <TGeoMatrix.h>
12 #include <TGeoManager.h>
13 #include <TGeoPhysicalNode.h>
14 #include "AliFMDGeometry.h"
16 //____________________________________________________________________
18 AliFMDSurveyToAlignObjs::GetUnitFactor() const
20 // Returns the conversion factor from the measured values to
22 if (!fSurveyObj) return 0;
23 TString units(fSurveyObj->GetUnits());
24 if (units.CompareTo("mm", TString::kIgnoreCase) == 0) return .1;
25 else if (units.CompareTo("cm", TString::kIgnoreCase) == 0) return 1.;
26 else if (units.CompareTo("m", TString::kIgnoreCase) == 0) return 100.;
30 //____________________________________________________________________
32 AliFMDSurveyToAlignObjs::GetPoint(const char* name,
34 TVector3& error) const
36 // Get named point. On return, point will contain the point
37 // coordinates in centimeters, and error will contain the
38 // meassurement errors in centimeters too. If no point is found,
39 // returns false, otherwise true.
40 if (!fSurveyPoints) return kFALSE;
42 Double_t unit = GetUnitFactor();
43 if (unit == 0) return kFALSE;
45 TObject* obj = fSurveyPoints->FindObject(name);
46 if (!obj) return kFALSE;
48 AliSurveyPoint* p = static_cast<AliSurveyPoint*>(obj);
49 point.SetXYZ(unit * p->GetX(),
52 error.SetXYZ(unit * p->GetPrecisionX(),
53 unit * p->GetPrecisionY(),
54 unit * p->GetPrecisionZ());
58 //____________________________________________________________________
60 AliFMDSurveyToAlignObjs::CalculatePlane(const TVector3& a,
67 // Vector a->b, b->c, and normal to plane defined by these two
69 TVector3 ab(b-a), bc(c-b);
71 // Normal vector to the plane of the fiducial marks obtained
72 // as cross product of the two vectors on the plane d0^d1
73 TVector3 nn(ab.Cross(bc));
74 if (nn.Mag() < 1e-8) {
75 Info("CalculatePlane", "Normal vector is null vector");
79 // We express the plane in Hessian normal form.
83 // where n is the normalised normal vector given by
85 // n_x = a / l, n_y = b / l, n_z = c / l, p = d / l
87 // with l = sqrt(a^2+b^2+c^2) and a, b, c, and d are from the
88 // normal plane equation
90 // ax + by + cz + d = 0
93 TVector3 n(nn.Unit());
94 // Double_t p = - (n * a);
96 // The center of the square with the fiducial marks as the
97 // corners. The mid-point of one diagonal - md. Used to get the
98 // center of the surveyd box.
102 // The center of the box.
103 TVector3 orig(md - depth * n);
108 // Normalize the spanning vectors
109 TVector3 uab(ab.Unit());
110 TVector3 ubc(bc.Unit());
112 for (size_t i = 0; i < 3; i++) {
113 rot[i * 3 + 0] = ubc[i];
114 rot[i * 3 + 1] = uab[i];
115 rot[i * 3 + 2] = n[i];
120 //____________________________________________________________________
122 AliFMDSurveyToAlignObjs::FitPlane(const TObjArray& points,
123 const TObjArray& errors,
124 Double_t /* depth */,
128 Int_t nPoints = points.GetEntries();
130 AliError(Form("Cannot fit a plane equation to less than 4 survey points, "
131 "got only %d", nPoints));
136 // Loop and fill graph
137 for (int i = 0; i < nPoints; i++) {
138 TVector3* p = static_cast<TVector3*>(points.At(i));
139 TVector3* e = static_cast<TVector3*>(errors.At(i));
141 if (!p || !e) continue;
143 g.SetPoint(i, p->X(), p->Y(), p->Z());
144 g.SetPointError(i, e->X(), e->Y(), e->Z());
147 // Check that we have enough points
149 AliError(Form("Only got %d survey points - no good for plane fit",
154 // Next, declare fitting function and fit to graph.
155 // Fit to the plane equation:
157 // ax + by + cz + d = 0
161 // z = - ax/c - by/c - d/c
163 TF2 f("plane", "-[0]*x-[1]*y-[2]",
164 g.GetXmin(), g.GetXmax(), g.GetYmin(), g.GetYmax());
167 // Now, extract the normal and offset
168 TVector3 nv(f.GetParameter(0), f.GetParameter(1), 1);
169 TVector3 n(nv.Unit());
170 Double_t p = -f.GetParameter(2);
172 // Create two vectors spanning the plane
173 TVector3 a(1, 0, f.Eval(1, 0)-p);
174 TVector3 b(0, -1, f.Eval(0, -1)-p);
175 TVector3 ua(a.Unit());
176 TVector3 ub(b.Unit());
177 // Double_t angAb = ua.Angle(ub);
178 // PrintVector("ua: ", ua);
179 // PrintVector("ub: ", ub);
180 // std::cout << "Angle: " << angAb * 180 / TMath::Pi() << std::endl;
182 for (size_t i = 0; i < 3; i++) {
183 rot[i * 3 + 0] = ua[i];
184 rot[i * 3 + 1] = ub[i];
185 rot[i * 3 + 2] = n[i];
188 // The intersection of the plane is given by (0, 0, -d/c)
197 //____________________________________________________________________
199 AliFMDSurveyToAlignObjs::MakeDelta(const char* path,
202 TGeoHMatrix& delta) const
204 // Make delta transformation
205 if (!gGeoManager) return kFALSE;
206 if (!gGeoManager->cd(path)) return kFALSE;
209 TGeoMatrix* global = gGeoManager->GetCurrentMatrix();
211 PrintRotation(Form("%s rot:", global->GetName()),global->GetRotationMatrix());
212 PrintVector(Form("%s trans:", global->GetName()),global->GetTranslation());
215 return MakeDelta(global, rot, trans, delta);
218 //____________________________________________________________________
220 AliFMDSurveyToAlignObjs::MakeDelta(TGeoMatrix* global,
223 TGeoHMatrix& delta) const
225 // Make delta transformation
226 TGeoHMatrix* geoM = new TGeoHMatrix;
227 geoM->SetTranslation(trans);
228 geoM->SetRotation(rot);
230 delta = global->Inverse();
231 delta.MultiplyLeft(geoM);
236 //____________________________________________________________________
238 AliFMDSurveyToAlignObjs::GetFMD1Plane(Double_t* rot, Double_t* trans) const
241 // The possile survey points
242 TVector3 icb, ict, ocb, oct, dummy;
244 if (!GetPoint("V0L_ICB", icb, dummy)) missing++;
245 if (!GetPoint("V0L_ICT", ict, dummy)) missing++;
246 if (!GetPoint("V0L_OCB", ocb, dummy)) missing++;
247 if (!GetPoint("V0L_OCT", oct, dummy)) missing++;
249 // Check that we have enough points
251 AliWarning(Form("Only got %d survey points - no good for FMD1 plane",
257 const char* lidN = "FMD1_lid_mat0";
258 TGeoMatrix* lidM = static_cast<TGeoMatrix*>(gGeoManager->GetListOfMatrices()
261 AliError(Form("Couldn't find FMD1 lid transformation %s", lidN));
265 const Double_t* lidT = lidM->GetTranslation();
266 Double_t lidZ = lidT[2];
267 Double_t off = lidZ-3.3;
272 if (!CalculatePlane(ocb, icb, ict, off, trans, rot)) return kFALSE;
273 // PrintRotation("FMD1 rotation:", rot);
274 // PrintVector("FMD1 translation:", trans);
279 //____________________________________________________________________
281 AliFMDSurveyToAlignObjs::DoFMD1()
284 Double_t rot[9], trans[3];
285 if (!GetFMD1Plane(rot, trans)) return kFALSE;
286 // const char* path = "/ALIC_1/F1MT_1/FMD1_lid_0";
288 // TGeoHMatrix delta;
289 TGeoTranslation global(0,0,324.670);
290 if (!MakeDelta(&global, rot, trans, fFMD1Delta))
293 // PrintRotation("FMD1 delta rotation:", fFMD1Delta.GetRotationMatrix());
294 // PrintVector("FMD1 delta translation:", fFMD1Delta.GetTranslation());
299 //____________________________________________________________________
301 AliFMDSurveyToAlignObjs::GetFMD2Plane(Double_t* rot, Double_t* trans) const
304 // The possible survey points
305 const char* names[] = { "FMD2_ITOP", "FMD2_OTOP",
306 "FMD2_IBOTM", "FMD2_OBOTM",
307 "FMD2_IBOT", "FMD2_OBOT",
309 const char** name = names;
314 // Loop and fill graph
317 if (!GetPoint(*name++, p, e)) continue;
319 points.Add(new TVector3(p));
320 errors.Add(new TVector3(e));
322 if (points.GetEntries() < 4) {
323 AliWarning(Form("Only got %d survey points - no good for FMD2 plane",
324 points.GetEntries()));
328 return FitPlane(points, errors, 0, trans, rot);
331 #define M(I,J) rot[(J-1) * 3 + (I-1)]
332 //____________________________________________________________________
334 AliFMDSurveyToAlignObjs::DoFMD2()
337 Double_t rot[9], trans[3];
338 if (!GetFMD2Plane(rot, trans)) return kFALSE;
339 // PrintRotation("FMD2 rotation:", rot);
340 // PrintVector("FMD2 translation:", trans);
342 for (int i = 0; i < 3; i++) {
343 for (int j = 0; j < 3; j++) {
344 rot[i*3+j] = (i == j ? 1 : 0);
347 trans[0] = trans[1] = 0;
349 // PrintRotation("FMD2 rotation:", rot);
350 // PrintVector("FMD2 translation:", trans);
352 // TGeoHMatrix delta;
353 if (!MakeDelta("/ALIC_1/F2MT_2/FMD2_support_0/FMD2_back_cover_2",
354 rot, trans, fFMD2Delta)) return kFALSE;
356 // PrintRotation("FMD2 delta rotation:", fFMD2Delta.GetRotationMatrix());
357 // PrintVector("FMD2 delta translation:", fFMD2Delta.GetTranslation());
362 //____________________________________________________________________
364 AliFMDSurveyToAlignObjs::Run()
366 AliFMDGeometry* geom = AliFMDGeometry::Instance();
368 geom->InitTransformations();
374 //____________________________________________________________________
376 AliFMDSurveyToAlignObjs::CreateAlignObjs()
378 TClonesArray& array = *fAlignObjArray;
379 Int_t n = array.GetEntriesFast();
381 if (!fFMD1Delta.IsIdentity()) {
382 new (array[n++]) AliAlignObjParams("FMD1/FMD1_T", 0, fFMD1Delta, kTRUE);
383 new (array[n++]) AliAlignObjParams("FMD1/FMD1_B", 0, fFMD1Delta, kTRUE);
385 if (!fFMD2Delta.IsIdentity()) {
386 new (array[n++]) AliAlignObjParams("FMD2/FMD2_T", 0, fFMD2Delta, kTRUE);
387 new (array[n++]) AliAlignObjParams("FMD2/FMD2_B", 0, fFMD2Delta, kTRUE);
394 //____________________________________________________________________
396 AliFMDSurveyToAlignObjs::PrintVector(const char* text, const TVector3& v)
399 Double_t va[] = { v.X(), v.Y(), v.Z() };
400 PrintVector(text, va);
402 //____________________________________________________________________
404 AliFMDSurveyToAlignObjs::PrintVector(const char* text, const Double_t* v)
408 << std::setw(15) << v[0]
409 << std::setw(15) << v[1]
410 << std::setw(15) << v[2]
415 //____________________________________________________________________
417 AliFMDSurveyToAlignObjs::PrintRotation(const char* text, const Double_t* rot)
419 // Print a rotation matrix
420 std::cout << text << std::endl;
421 for (size_t i = 0; i < 3; i++) {
422 for (size_t j = 0; j < 3; j++)
423 std::cout << std::setw(15) << rot[i * 3 + j];
424 std::cout << std::endl;
428 //____________________________________________________________________